Residual sum of squares in a regression I understand that in a linear regression model, the residual sum of squares will either remain same or fall with the addition of a new variable. 
What if the two models were 
$$
  I \colon y_i=\beta_0 + \beta_1 x_{1i}+\epsilon_{i}
$$
and
$$
  II \colon  y_i = \beta_0 + \beta_1 x_{1i} + \beta_2 {x_{1i}}^{2} + \epsilon_i
$$
Then, will the residual sum of squares of model 2 be less than or equal to that of model 1?
 A: While it is true that your second model is technically also a linear regression model, I don't think that's the salient point here. 
The point here is that $M_2(x) = M_1(x) + f(x|\theta)$
in which $M_{1}(x)$ is "Model 1" and $M_{2}(x)$ is "Model 2".
Provided there exists a $\theta$ such that $f(x|\theta)=0 \forall x$ (see comments for a discussion of this),
it must be the case, that on your training data, the sum of residuals for model 2 is less than or equal to your residuals for model 1.
Why is this the case? Imagine you wish to train Model 2, but you initialise it such that $M_{2}(x)=M_{1}(x)$. In this case, that means you choose $\beta_{0}$ and $\beta_{1}$ to take the values they take for $M_{1}(x)$, and you choose $\theta$ s.t. $f(x|\theta)=0 \forall x$. You then start gradient descent. One of two things will happen. Either, gradient descent will decide that there is a combination of parameters which decreases the loss (your sum of residuals) further, or it will determine you're in a local minimum and exit (this is the case for which the sums of residuals are the same for both models). Remember that just because your training loss for your second model is less than or equal to that for model 1, doesn't mean it's a better model, you could be overfitting and your test loss could be worse for model 2. 
A: One is not supposed to believe this the first she sees it, but that model II you posted is a linear regression.
The linearity has nothing to do with the variables. “Linear” refers to linearity in the parameters. Now if you were to take $(\beta_2x_{1i})^2$, then you would have a nonlinear transformation of a parameter and a nonlinear regression. As it stands, though, squaring a variable doesn’t change your regression to nonlinear.
For a nonlinear regression, the answer is that in-sample performance can get worse when you add another parameter. The idea in linear regression is that, if all else fails, the new parameter estimate can be 0 and the old parameter estimates the same as before, so performance will be identical. 
In the nonlinear case, consider a model (might be linear, might not be) that has perfect performance. For whatever reason, you decide to add an $e^{\beta x}$ term. This is a nonlinear transformation of the $\beta$ parameter, so the regression is now guaranteed to be nonlinear. However, the entire $e^{\beta x}$ term is always positive, no matter the value of $\beta$. Therefore, the fit will be different than the perfect model, so performance will drop.
(This argument isn’t perfect, since the new model with the exponential term might also happen to fit the data perfectly, but I think this illustrates the point that setting the new parameter to zero need not matter when there is a nonlinear transformation of a parameter.)
Edit
This only applies to in-sample data. Out-of-sample, even if the higher-parameter model does better in-sample, it may just have fit to coincidences in the data, not the real trend, meaning that you may (will) do poorly when you make predictions about additional values.
Edit 2 (skip down to #4)
I'm not convinced that this post is wrong, but I have a different argument. Unlike linear regression, adding a parameter to a nonlinear regression does not have to involve adding a term. 
Indeed, $y = 2^{\beta x}$ to $y=\dfrac{2^{\beta x}}{\beta}$ is something like adding a parameter, even if that parameter happens to be the same parameter.
However, say that $y = 2^{2 x}$ is a perfect fit, such as when  $x = (1,2,3,4,5)$ and $y=(4,16,64,256,1024)$. 
Then $y = \dfrac{2^{2 x}}{2}$ will not be a perfect fit. 
However, $y = \dfrac{2^{1x}}{1}$ is not a perfect fit either. No matter what we do, we either have the wrong exponent or the wrong denominator.
Edit 3
Working through these, I'm not sure that I've given an example. However, there are so many funky nonlinear functions that I am confident that there is one. Maybe it would be something like $e^{1+\beta^2} * I(x_i = 1)$ where $I(x_i = 1)$ is an indicator function that $x_i = 1$.
Edit 4
Let $X=(1,2,3,4,5,6,7,8,9,10)$ and $Y=(2,3,4,5,6,7,8,9,10,11)$.
Then the perfect fit of $y_i = \beta_0 + \beta_1x_i$ has $\beta_0 = 1, \beta_1=1$.
Now consider the fit of $y_i = \beta_0 + \beta_1x_i + \dfrac{\beta_1^{\beta_2}}{1-\beta_1}$.
We need $\beta_1=1$ to get the right slope, since nothing else affects the slope. However, then we divide by zero, so we can't ever have the right slope. Therefore, we won't ever get a perfect fit.
Now we have nested models (as whuber mentions is important for doing this sort of comparison, except using $\beta_1$ twice means that it actually isn't nested), but the more complicated model will have a worse fit even in sample.
(There should be hats (e.g. $\hat{\beta}$) all over the place in this post, yes.)
A: Linear models
We first assume that the focus of the question is the two models posted. In those two models (which are actually linear in the coefficients, not non-linear) we can take a constructive approach noting due to the nesting of the models that the residual sum of squares (RSS) of model 1 minus the RSS of model 2 equals the squared length of the projection of the residual vector of model 1 on the range space of the X matrix of model 2.
Since a squared value cannot be negative RSS for model 1 must greater or equal to RSS for model 2.
Nonlinear models
Now we consider nested nonlinear models. Since linear models are a subset of nonlinear models this applies to them too. We require that model 1 is model 2 subject to an additional restriction on the latter's coefficients.  Realize that the RSS of model 2 is the minimum value of $ || y - f(X,  \beta) ||^2 $ over $\beta$ for some fixed known function $f$.  For example, if the model is linear then $f(X, \beta) = X \beta$.  The RSS of model 1 is the minimum of that same objective function but with an additional constraint.  e.g. $\beta_2 = 0$. Since model 1 has an additional constraint model 1's objective after optimizing can be no lower than that of model 2.
Example - linear models
For the linear models discussed we can verify the above with a particular example using R code:
# test inputs
set.seed(123)
x <- 1:5
y <- 1 + x + x*x + rnorm(5)

# compute the two models
fm1 <- lm(y ~ x)
fm2 <- lm(y ~ x + I(x^2))

# projection of residual vector of model 1 on range of X matrix of model 2
# Use fact that fitted values of a regression is the projection of the
# response vector on the range of the X matrix.
p <- fitted(update(fm2, resid(fm1) ~ .))

# squared length of p
crossprod(p)
##          [,1]
## [1,] 7.402135

# difference between RSS values is same
deviance(fm1) - deviance(fm2)
## [1] 7.402135

Also, if $X$ is the model matrix of the second model and $X_0$, $X_1$ and $X_2$ are its three columns and if the first two columns of $X$ form the model matrix of the first model and if $X_2$ has unit length and is orthogonal to the other two columns then the squared length of the projection discussed above can be expressed as $(X_2'y)^2$ so:
fm1a <- lm(y ~ poly(x, 1))
fm2a <- lm(y ~ poly(x, 2))

# last row of cross prod matrix shows 3rd col has
# length 1 and is orthogonal to others
X <- model.matrix(fm2a)
crossprod(X)[3, ]
##   (Intercept)   poly(x, 2)1   poly(x, 2)2 
##  1.110223e-16 -1.110223e-16  1.000000e+00 

crossprod(X[, 3], y)^2
##          [,1]
## [1,] 7.402135

# difference between RSS values is same
deviance(fm1a) - deviance(fm2a)
## [1] 7.402135

