How to motivate the definition of $R^2$ in `sklearn.metrics.r2_score`? TLDR: What motivates the definition of $R^2$ in the Python function sklearn.metrics.r2_score?
DETAILS
The Python machine learning package sklearn implements an $R^2$ using the following formula.
$$
R^2=1-\left(\dfrac{
\overset{N}{\underset{i=1}{\sum}}\left(
y_i-\hat y_i
\right)^2
}{
\overset{N}{\underset{i=1}{\sum}}\left(
y_i-\bar y
\right)^2
}\right)
$$
This is fine for in-sample assessments of $R^2$, but when it comes to out-of-sample assessments, the definition seems unmotivated and lacking in meaning.
In simple linear regression (in-sample), there are multiple equivalent definitions of $R^2$.

*

*Squared correlation between the feature and the outcome

*Squared correlation between the outcome and the predictions

*Proportion of variance explained

*Comparison of the square loss incurred by the model to the square loss incurred by a baseline model

(There may be more, and perhaps it is me only thinking of these four that prevents me from seeing what this Python function means.)
In regression models with multiple features, definition 1 does not make sense, so it is not a contender for a generalization of $R^2$ to out-of-sample assessments in complicated situations that are likely to have multiple features.
Out-of-sample, squaring the correlation misses prediction bias. For instance, $y=(1,2,3)$ is perfectly correlated with $\hat y=(11,12,13)$, yet $\hat y$ is composed of terrible predictions of $y$. While this kind of bias cannot happen in-sample in ordinary least squares linear regression, out-of-sample, all bets are off (and that’s without even getting into what could go wrong in fancier regression models like support vectors or neural networks). Consequently, definition 2 is out as a generalization of $R^2$ to out-of-sample assessments.
Except in special situations, the usual methods for calculating $R^2$ do not correspond to the proportion of variance explained, and if we twist the definition to force it, we usually lose the connection to the square loss that we aim to minimize. Thus, definition 3 does not seem like a good contender for a generalization of $R^2$.
Finally, definition 4 makes sense. We have some kind of baseline model (naïvely predict $\bar y$ every time, always using the marginal mean as our guess of the conditional mean) and compare our predictions to the predictions made by that baseline model. To draw an analogy to flipping a coin, if someone guesses which side will land up and gets correct predictions less than half the time, that person is a poor predictor. If they right more than half the time, they are at least improving somewhat upon the naïve “Gee, I don’t know how it’ll land, so I guess I’ll just say heads every time (or tails, or alternate between the two) and get it right about half the time.”
(Put bluntly, why pay a data scientist or statistician a lot of money to predict conditional means when you could do better just by predicting one number every time that you can calculate by AVERAGE(A:A) in Excel?)
In-sample, I am totally on board with the Python function. Out-of-sample, I have a problem. In the $R^2$ formula above, the Python implementation uses the $\bar y$ from the given data. That is, if you input an out-of-sample $y$ and the corresponding predictions $\hat y$, the function uses the out-of-sample $\bar y$ in the denominator.
This corresponds to comparing your model to a model that predicts the out-of-sample mean every time. However, we cannot have access to such a model, since it would require knowledge of the true out-of-sample values.
QUESTION: What motivates this definition that is used in sklearn?
An argument about ease of use does not seem legitimate to me (even if that is why the function is defined this way), because the following is even easier to implement.
def r2_score():
    return(0.5)
    # or return(np.random.uniform(0, 1, 1))

While this is an extreme example, the function does not return useful information about the regression model. Consequently, making a function easy to use at the expense of returning unhelpful information seems like a poor design or a function that we should not use.



*

*Squared correlation between the feature and the outcome


That would be the case if you have a single feature and the model is linear regression.



*Squared correlation between the outcome and the predictions


Same as above, but it will hold also if there are more features.



*Proportion of variance explained


It tells us the proportion of the variance explained, but only for the linear regression.



*Comparison of the square loss incurred by the model to the square loss incurred


Again, for linear regression using the formula from Scikit-learn is equivalent to the as we can decompose the squared error to TSS = ESS + RSS and get the equivalent formulation.
So in a sense, all the formulations are the same, just have varying degrees of generality.
As for calculating $R^2$ on the test set, you can check this thread and if you search through CrossValidated.com and Scikit-learn GitHub issues, you can find many discussions considering the Scikit-learn's choice of test set mean in the denominator as controversial. As you can learn from this discussion, one of the problems with using the train set mean for the test set $R^2$ is the API, where the metrics usually are defined as metric(y_true, y_pred) and the same interface can be used regardless if it is training or test set. It would also mean that during evaluation time you would need to have access to the training data, which may not be possible in some setups where there is a hard split of the data between training time and validation time.
Notice also that $R^2$ comes from statistics, where we are usually interested in in-sample metrics, so the derivations would as well regard the train set $R^2$.
