Why can image denoising be expressed as an energy minimization problem?

I am interested in this paper by (Ulyanov et al., 2017), but I don't have much knowledge about image denoising. I don't understand why image denoising can be expressed as an energy minimization problem. Where does the following formula come from and what it does it mean?

To demonstrate the power of this parametrization, we consider inverse tasks such as denoising, super-resolution and inpainting. These can be expressed as energy minimization problems of the type $$x^∗ = \min_x E(x; x_0) + R(x), \text{ (1)}$$

where $E(x; x_0)$ is a task-dependent data term, $x_0$ the noisy/low-resolution/occluded image, and $R(x)$ a regularizer.

Ulyanov, Dmitry, Andrea Vedaldi, and Victor Lempitsky. "Deep Image Prior." arXiv preprint arXiv:1711.10925 (2017).

The very short answer is that energy minimisation is a general, unifying framework for both probabilistic and deterministic machine learning models. In this very general setting, energy does not have much to do with a physical notion of energy and it is essentially a probability which has not been normalised to fall within $[0, 1]$.
In image denoising there is a bit of a deeper connection with physics which partly explains where this energy minimisation language comes from. Let's only consider binary images $X$ of size $n \times n$, so $X \in \{-1,1\}^N$. Suppose that we believe that the value of a pixel $s$ should depend on it's neighbours, that is the pixels which share an edge with $s$. Let $X^+_s$ and $X^−_s$ denote the number of neighbours of $s$ that take positive and negative values, then we can formalise our assumption of spatial dependency as follows $$P(X_s = +1|X_t, t\text{ is a neighbour of s}) = \frac{exp(2 \beta (X^+_s − X^−_s))}{1 + exp(2\beta(X^+_s − X^−_s))}$$ where $\beta > 0$. The above model also occurs in statistical physics where it is called the Ising model of ferromagnetism. The parameter $\beta$ is the inverse temperature which is very closely related to energy.
So if the noise is spatially independent, then computing the most likely configuration $X^*$ that satisfies the above relation will effectively denoise our image. Computing the most likely configuration is essential done by minimising the energy of the Ising model.
The equation $x^∗ = \min_x E(x; x_0) + R(x)$ is a pretty generalised formulation of the denoising problem. I'll give an example of a very common method which fits into this framework. Suppose we have an image $X$ which has been corrupted with additive Gaussian noise. That is we have $$Y = X + \epsilon, \text{ where } \epsilon \sim N(0, \sigma)$$ where $Y$ is the observed noisy image and $\epsilon$ is the noise. A popular way to remove the noise $\epsilon$ is with an algorithm called basis pursuit. To denoise with basis pursuit we choose a prior that $X$ is generate from a sparse, redundant set of visual patterns $D$ which we are given. Then to recover our original image, we split $Y$ up to into a series of patches $Y_{ij}$ and try to minimise $$\min_\alpha \{ ||Y_{ij} - D \alpha||^2_2 - \lambda||\alpha||_1 \}$$ for each patch. Here $E(x; x_0) = ||Y_{ij} - D \alpha||^2_2$ and $R(x)= \lambda||\alpha||_1$. We could also choose different prior about what generated $X$ and this would result in a different $R(x)$.