Generating nonlinearities in renormalization group In renormalization group (RG) calculations as performed in statistical physics (for example for Landau-Ginzburg theory - often a la Wilson), the first step is to coarse-grain the theory by integrating out the fast degrees of freedom. Starting from some specified interactions, performing this step in practice often leads to new forms of nonlinear terms in the theory‌. It is then argued for a consistent calculation, the starting Hamiltonian should contain all possible interaction terms that are allowed by the symmetries and constraints of the problem at hand.
Although this is a very standard approach, it leaves a question: exactly how new interaction terms are "generated" by coarse-graining a given set of interactions?
A backward argument could be made on the basis of symmetries for when new terms will not appear in the RG ("if one/any interaction is not going to be generated at any level of perturbation, there must be a symmetry preventing their generation");  I'm however mostly looking for an intuitive way of understanding the generation of new nonlinearties from old ones, rather than lack thereof.
 A: To start, note that the generation of new interactions has everything to do with coarse-graining, it only appears in RG calculations because coarse-graining is the first step of an RG procedure. The essential issue is that if one has a system comprising multiple interacting random variables and you want to describe the behavior of just a subset of those random variables, then that description of the subset must contain effective interactions between random variables that do not necessarily have any direct interaction. This is true whether you have a field theory, multivariate distribution, or even a bivariate distribution. I'll leave aside the issue of relevant and irrelevant interactions here, which come up in the second step of the RG procedure, rescaling.
I always like to go back to simpler models in statistical mechanics to think about these things. As a simple example, let's consider the classical Ising model and a one-step coarse graining in real space by integrating out every other spin. In the full Hamiltonian for the model in which all of the spins are observed, suppose each spin only interacts with its nearest neighbors. e.g., in $d=1$
$$\mathcal H(\mathbf{s}) = -J \sum_{i} s_i s_{i+1} - h s_i$$
with $P(\mathbf{s}) = \exp(-\beta \mathcal H(\mathbf{s}))/Z$ and $Z$ the normalization factor. The coarse-grained Hamiltonian for the 1d example would be defined, say, by summing over all of the even sites:
$$\exp(-\beta\mathcal H'(\mathbf{s})) = \sum_{s_{2n}} \exp(-\beta \mathcal H(\mathbf{s})).$$
Note that if we know the full Hamiltonian, we can in principle calculate any statistical moment of the distribution: mean orientation of each spin, spin-spin correlations, etc. Our coarse-grained Hamiltonian must give the exact same values for these statistics; e.g., $\langle s_1 s_5 \rangle$ must be the same value whether we compute it using the full Hamiltonian $\mathcal H(s)$ or the coarse-grained Hamiltonian $\mathcal H'(s)$.
We know it cannot possibly be the case that our coarse-grained Hamiltonian contains only independent spins, since then all of the higher order statistics like correlations would be $0$. It must be the case that the coarse-grained Hamiltonian has developed effective interactions between spins that did not directly interact in the fully resolved model. One might hope that these effective interactions only couple what were formerly next-nearest neighbor spins, but often that is not exactly the case: any interaction allowed by symmetry has developed, although the relative magnitudes of these interactions are not all the same---far away spins would have weaker interactions than nearby spins in this one-step coarse-grained model.
As an even simpler example, let's consider a zero-dimensional model containing just two random variables $x$ and $y$, both taking values on $\mathbb{R}$. Suppose the probability density describing the interactions between these two random variables is
$$\rho(x,y) \propto \exp\left(-x^2 - x^4 - y^4 - x^2y^2\right).$$
Some things to note about this distribution: it is even in both $x$ and $y$, and it does converge. Suppose we only care about the variable $x$ (say because it is easier to observe), so we would be just as happy with a model for only $x$. The proper way to obtain this model is to marginalize out $y$. In this case one can perform the quartic integral in terms of an infinite series:
$$\int_{-\infty}^\infty dy~e^{-x^2y^2 - y^4} = \frac{1}{2}\sum_{p=0}^\infty \frac{(-x^2)^{p}}{p!} \Gamma\left(\frac{2p+1}{4}\right) = e^{\frac{x^4}{8}} |x| K_{1/4}\left(\frac{x^4}{8} \right);$$
for completeness I've added the representation of the series in terms of the modified Bessel function $K_\alpha(x)$, but the series is more useful for our general purpose. Now let us write the marginalized distribution in the form
$$\rho(x) \propto \exp\left(-x^2 - x^4 + \ln\left[\frac{1}{2}\sum_{p=0}^\infty \frac{(-x^2)^{p}}{p!} \Gamma\left(\frac{2p+1}{4}\right) \right] \right).$$
Thus, we see that even in this $0d$ example marginalizing out one of the random variables has resulted in a marginal distribution for $x$ that has not only modified the coefficients of the existing terms (if we extract the $x^2$ and $x^4$ terms by expanding the logarithm in a formal series), but has also "generated" even powers of all higher orders---note that there are no odd powers, consistent with the symmetry of the full distribution. These higher order interactions must be there to ensure that the statistical moments $\langle x^m \rangle = \int dx~\rho(x) x^m $ match the moments obtained by calculating $\langle x^m \rangle = \int dx dy~\rho(x,y) x^m$. (You can check numerically that, e.g., the second moment evaluates to $\sim 0.21$ using both the full and marginalized distributions. For the latter, the representation of the series in terms of the Bessel function is useful in, e.g., Mathematica).
A: A simple way to see intuitively how different coupling may be generated as a result of the coarse-graining is the following example, adapted from Maris and Kadanoff pedagogical introduction to RG (Maris, H. J., & Kadanoff, L. P. (1978). Teaching the renormalization group. American Journal of Physics, 46(6), 652-657).
Let's use the Ising model. A convenient coarse-graining can be obtained by a partial trace where one degree of freedom every two is summed.
In 1D, this is equivalent to eliminating the b, d, and sites in the following figure. The remaining a', c', e',... sites interact as an effect of the previous level interactions between a and, band c and so on. The original first neighbor interaction is transformed into a new first neighbor interaction, although with a different coupling constant, due to the summing procedure.

In 2D things are qualitatively different.
A similar elimination corresponds to summing on the small circles (sites) in the following picture.

Summing the direct interactions HC, HG, HM, HI induces additional interactions between the survival sites C, G, M, I, according to the following figure:

The new interactions correspond to nearest-neighbor interactions in the new lattice (red lines), next-nearest-neighbor interactions (green lines), and an irreducible 4-spin interaction among the C, G, M, I sites.
A: The rough idea is this: Say you know that a particular theory describes real-world observations for all energies $E \in [0, \Lambda_0]$. For the sake of argument let this particular theory be a scalar QFT with only the $\phi^4$ interaction term. This means that $\phi(x)$ is defined only for $E \in [0, \Lambda_0]$.
What happens when you integrate out the high energy Fourier modes of $\phi(x)$ in the range $E \in[\Lambda, \Lambda_0]$, where $\Lambda < \Lambda_0$?
The new theory (with fields $\tilde \phi(x)$) that you get from this procedure now describes physics only for energy scales $E \in [0, \Lambda]$. But since we got this new theory from the old one, it must still reproduce all the predictions that the previous theory defined for $E \in [0, \Lambda_0]$ did.
The only way this is possible and consistent is if the new theory has new interaction terms and different coupling constants. And, as you point out, these new terms can be of a general form which the symmetries allow.
Typically we say that the new theory is an effective field theory describing the old physics. The renormalisation group precisely tells you which operators will appear in the effective theory describing the (old) high-energy physics. TL;DR new terms need to be introduced because the coarse-grained theory must consistently describe the original theory.
