Does $R^2$ increase with the number of variables no matter what model we are adjusting? I understand this is true for OLS, but i am not sure if this is true for every other model, like ridge or a NN.
 A: There are less common estimators that we can use and might want to use under some circumstances. For instance, if we estimate the parameters of a linear model by minimizing something other than square loss (sum of squared residuals), we might wind up with the square loss increasing as we add parameters, since the optimization did not know to try to minimize square loss.
In the simulation below, I have an estimator that minimizes a quantile loss to estimate the coefficients. The result is the smaller model has a smaller sum of squared residuals and, consequently, a larger $R^2$.
library(quantreg)
set.seed(2022)
N <- 1000
x1 <- rnorm(N)
x2 <- rnorm(N)
yh <- x1
e <- rnorm(N, 0, 1)
y <- yh + e
L <- lm(y ~ x1 + x2)
X1 <- cbind(x1)
X2 <- cbind(x1, x2)

quantile_r2 <- function(X, y){
  
  L <- quantreg::rq(y ~ X, tau = 0.6)
  
  preds <- predict(L)
  
  SSTotal <- sum(
    (
      y - mean(y)
    )^2
  )
  
  SSResiduals <- sum(
    (
      y - preds
    )^2
  )
  
  R2 <- 1 - SSResiduals/SSTotal
  
  return(R2)
}
quantile_r2(X1, y) - quantile_r2(X2, y)

$$
R^2 = 1-\dfrac{
\underset{i=1}{\overset{n}{\sum}}\left(
y_i - \hat y_i
\right)^2
}{
\underset{i=1}{\overset{n}{\sum}}\left(
y_i - \bar y
\right)^2
}
$$
Overall, it depends on what exactly you're doing if you can claim that more parameters necessarily results in higher in-sample $R^2$.
A: The question is a bit vaguely formulated (how is one increasing the number of variables?, what sort of loss function is your model trying to minimize?). In the following result, I try to state a set of precise conditions under which the answer to your question is "yes". I then discuss why my result is true under the stated conditions as well as when my stated conditions fail to give you a sense for what sort of "pathologies" as alluded to in @whuber's comment one would need for $R^2$ to not increase.
Claim: Let $X_1$ be a vector of predictors and let $X_2$ be a vector of predictors containing $X_1$ as a sub-vector. Let $Y$ be the dependent variable whose value is being predicted. Let $\theta_1\in \mathcal H_1$ be the parameters of a model using $X_1$ as predictors such that when the model parameter is $\theta_1$, our predictions of $Y$ given $X_1$ are $m_1(X_1;\theta_1)$ and let $\theta_2 \in \mathcal H = \mathcal H_1 \cup \mathcal H_2$ be the parameters of a model using $X_2$ as predictors such that when the model parameter is $\theta_2$, our prediction of $Y$ given $X_2$ is $m_2(X_2,\theta_2)$. Finally, assume that our models are nested in the sense that if $\theta_2\in\mathcal H_1$, then $m_1(X_1,\theta_2) = m_2(X_2,\theta_2)$ whenever $X_1$ has the same values as $X_2$ in all shared commponents. Finally, suppose that our "best" predictors are chosen via least squares, i.e. for $i=1,2$,
$$\theta_i^* = \mathrm{argmin}_{\theta_i}\, \mathbb E[(Y-m_i(X_i,\theta_i))^2]$$
Define the population $R^2$ of model $i$ by the quantity
$$R^2_i \equiv 1 - \frac{\mathbb E[(Y-m_i(X_i,\theta_i^*))^2]}{\mathrm{Var}(Y)}$$
Then $R_1^2 < R_2^2$.
The claim required a lot of setup to state, but the proof is simple. Note that when one is fitting a least-squares model, the objective being minimized is exactly the numerator of the object being subtracted in $R^2$, and the denominator depends only on properties of $Y$. Moreover, when fitting using predictors $X_2$, our assumption that the models are nested implies that a valid possible value of $\theta_2$ is one that effectively ignores the extra information in $X_2$ and predicts according to $\theta_1^*$. Thus, the fact that $\theta_2^*$ minimizes the objective implies that one can only do better with a richer set of predictors. Formally, by definition of minimization, we have
$$\mathbb E[(Y-m_2(X_2,\theta_2))^2] \leq \mathbb E[(Y-m_2(X_2,\theta_1^*))^2] = \mathbb E[(Y-m_1(X_1,\theta_1^*))^2]$$
and the rest follows from some simple algebra.
So now, what are the possible violations to my long list of assumptions. First, it is possible that $X_2$ has a larger number of variables than $X_1$, but there are some variables in $X_1$ that are highly predictive but not in $X_2$. In this case, model 2 might have a lower population $R^2$. Second, it is possible that the objective function used to choose the "best" $\theta_i$ is not the least-squares objective. This is essentially @Dave's example, and there again, there is no reason in general to assume that a $\theta_i$ chosen to fit some other objective happens to also fit the least-squares objective well too.
