Why does the rank order of models differ for R squared and RMSE? I am comparing $R^2$ and RMSE of different models. Interestingly, the rank ordering of the models with respect to $-R^2$ and RMSE is different and I do not understand why.
Here is an example in R:
library(caret) 

set.seed(0)
d<-SLC14_1(n=1000)

folds<-createMultiFolds(d$y,k=10,times=1)
tc<-trainControl(index=folds,returnResamp="all")
t1<-train(y~.,data=d,method="glmnet",trControl=tc) 
order(t1$results$RMSE)==order(-t1$results$Rsquared)

Output:
[1]  TRUE FALSE  TRUE FALSE FALSE FALSE FALSE FALSE  TRUE

Thus, the order if different for $-R^2$ suqared and $RMSE$.
The question is, why.
Let $SS_{res}$ be the sum of squared residuals $\sum(y_i-f_i)^2$. 
$RMSE$ is defined as $\sqrt{SS_{res}/n}$.
$R^2$ is defined as $1-SS_{res}/SS_{tot}$ where $SS_{tot}$ is $\sum(y_i-\overline{y})^2$. 
Since $SS_{res}=n*(RMSE)^2$, we can write $R^2$ as $1-n*(RMSE)^2/SS_{tot}$.
Since $n$ and $SS_{tot}$ are constant and the same for all models, $-R^2$ and $RMSE$ should strictly positively related. However, they are not since the ranking order is in practice not identical (see example code).
What is wrong with my argument?
 A: It's because caret calculates R-squared differently than you are.  See the answer to this question: How caret calculates R Squared.
To see it in your code, 
library(caret) 

set.seed(0)
d<-SLC14_1(n=1000)

folds<-createMultiFolds(d$y,k=10,times=1)
tc<-trainControl(index=folds,returnResamp="all",
             savePredictions = TRUE) # New option 
t1<-train(y~.,data=d,method="glmnet",trControl=tc) 
order(t1$results$RMSE)==order(-t1$results$Rsquared)

library(data.table)
preds <- data.table(t1$pred)
preds[, overall_mean := mean(obs), by = .(lambda, alpha, Resample)]

sum_sq <- preds[, .(SS_res = sum((obs - pred)^2),
                SS_tot = sum((obs - overall_mean)^2),
                n = .N,
                var = var(obs),
                Rsquared_corr = cor(obs, pred)^2),
            by = .(lambda, alpha, Resample)]
sum_sq <- sum_sq[, ':=' (RMSE_Julian = sqrt(SS_res / n),
                         Rsquared_Julian = 1 - (SS_res/SS_tot),
                         Rsquared_traditional = 1 - (SS_res/ ((n-1)*var) ))]
sum_sq <- merge(sum_sq, t1$resample, by = c("lambda", "alpha", "Resample"))
head(sum_sq)

Note the savePredictions = TRUE in the call to traincontrol().  In the final dataset, sum_sq, you can see your result, Rsquared_Julian matches Rsquared_traditional, but these don't match Rsquared_corr which does match the R-squared from caret, Rsquared.  
Also in your question, you assume n and SS_tot are constant, but that only holds true for a fold, not across all the cross-validations.
A: @grand_chat has the correct maths, I'm just growing in a comparative example to help illustrate what the issue is in different terms that will hopefully help understanding.
We're working with fractional terms here, similar to say miles per gallon. If we average mpg over set units of time we get very different results compared to over set units of fuel or distance. 
If we travel 10 minutes at 50 mph achieving 50mpg then 10 minutes at 60 achieving 30 mpg and we then want to calculate the average fuel efficiency for the journey.
Time based average (with one minute representing a unit of time) is $(50*10+30*10) /20 = 40 mpg$
But the distance we travel is $50/6 + 60/6 = 18.33 miles $ given that ten minutes is 1/6 th of an hour
The fuel we use is $(50/6)/50 +(60/6)/30=1/6+2/6 = 1/2 gallon$
This means our average mpg is in fact $18.33/(1/2)=36.66$
Because the total variance is different in every fold you would need to account for this in the averaging to maintain the monotonic relationship. Since it is present in the R2 calculation but not the RMSE then you can get rank switching by not accounting for the total variance in each fold 
