Sorta, yes, there is a classical analogue to renormalization groups that is native to classical field theory - at least in the case of scalar fields. I will demonstrate it here. Whether/how it can be dovetailed into the renormalization group framework in Quantum Field Theory, though, I will leave open, though it has one key point in common: renormalization pertains first and foremost to the constitutive coefficients in field theory. I will explain further.
Suppose we have a scalar field theory that is derived from an action principle associated with the action $S = \int 𝔏 d^nx$, where $n$ is the number of dimensions in the underlying geometry ($n = 4$ for space-time) and
$$d^n x = dx^0 ∧ dx^1 ∧ ⋯ ∧ dx^{n-1},$$
where $\left(x^0, x^1, ⋯, x^{n-1}\right)$ are the coordinates for the geometry. Assume the Lagrangian density $𝔏(q,v)$ depends only on a scalar field $q(x)$ and its gradients $v_μ = ∂_μq$. The coordinate differential operators will be denoted
$$∂_μ = \frac{∂}{∂x^μ},$$
as is mostly standard.
Finally, suppose that the geometry has a metric, given by the line element
$$g_{μν}dx^μdx^ν\quad \left(\text{i.e. }\sum_{0≤μ,ν<n} g_{μν}dx^μdx^ν\right),$$
where I will be using the summation convention, as shown here, below. For simplicity, we'll assume that this the metric coefficients are constant and that the metric is non-degenerate, with an inverse $g^{μν}$.
We want the Lagrangian to be invariant - at least with respect to transforms that preserve the metric. More generally, we'd want it to be invariant with respect to arbitrary coordinate transforms, but that would requiring adding the metric, itself (and its gradients) to the list of field components and to the Lagrangian. So, for simplicity, we'll just treat it as a constant.
If the Lagrangian is invariant, then it reduces to a function of the invariant combination of the field components. The following are the invariants:
$$I_0 = ½q^2,\quad I_1 = ½v·v = ½g^{μν}v_μv_ν.$$
Therefore, the Lagrangian density can be expressed as a function $𝔏\left(I_0, I_1\right)$ of the invariants.
Now... what does the Lagrangian density do? The field components $(q,v)$ play a role analogous to the Kinematic Variables in mechanics (i.e. positions, angles, and time rates of change thereof). Between them are Kinematic Relations, the role of which is played here by the field-gradient relation:
$$v_μ = ∂_μq.$$
The Lagrangian density generates the Dynamic Variables for the system being described by the action principle. In mechanics, this includes masses, momenta, angular momenta, energy, forces, torques and so on. They are generated as the derivatives with the Lagrangian density with respect to the field components. Here, the analogous relation is played by the following Lagrangian derivatives:
$$𝔉 = \frac{∂𝔏}{∂q},\quad 𝔓^μ = \frac{∂𝔏}{∂v_μ}.$$
The action principle produces, for them, a Dynamic Equation, given by the Euler-Lagrange equation, which here is:
$$∂_μ 𝔓^μ = 𝔉.$$
The form that the Kinematic and Dynamic equations take is independent of the Lagrangian density and does not convey anything empirical, in itself, other than that the system being described actually has these kinematic and dynamic variables for its description. All such systems, given in this way by an action principle for a Lagrangian density of the indicated type, must have these two equations.
Rather, the empirical content is contained in the form which the Lagrangian density takes and - by virtue of this - in the relations between the kinematic variables $(q,v)$ and dynamic variables $(𝔓,𝔉)$. These relations are the constitutive relations for the system in question - and this is precisely what the Lagrangian is there to specify.
The result of requiring that the Lagrangian density be invariant, and that it therefore be a function of the invariant combinations of field components, is that these relations can be completely encapsulated into a relatively small set of constitutive coefficients - namely: the derivatives of the Lagrangian density with respect to the invariants:
$$κ = \frac{∂𝔏}{∂I_0},\quad λ = \frac{∂𝔏}{∂I_1}.$$
Each of these coefficients, in general, is a function of the invariants:
$$κ = κ\left(I_0, I_1\right),\quad λ = λ\left(I_0, I_1\right),$$
such that the corresponding mixed-derivative relations for their derivatives holds:
$$\frac{∂^2𝔏}{∂I_1∂I_0} = \frac{∂^2𝔏}{∂I_1∂I_0}\quad⇒\quad\frac{∂κ}{∂I_1} = \frac{∂λ}{∂I_0}.$$
The encapsulation of the constitutive relations takes the following form, which can be derived by writing out the total differential for the Lagrangian density:
$$d𝔏 = κ dI_0 + λ dI_1 = κ q dq + λ g^{μν}v_ν dv_μ\quad⇒\quad \left(𝔉 = κ q,\quad 𝔓^μ = λ g^{μν} v_ν\right).$$
The forms that the constitutive relations take are also independent of the Lagrangian density! So, we can expand the previous remarks that the Kinematic and Dynamic equations are of a fixed Lagrangian-density-independent form to also include the constitutive relations. The only place where there's dependence is in the coefficients, themselves, and this is where the encapsulation takes place.
The resulting Euler Lagrange equation has the following fixed form:
$$∂_μ\left(λg^{μν}∂_νq\right) = κq,$$
or more briefly: $∂·(λ∂q) = κq$. This is a generalization of the Sturm-Lioville Problem in that (1) the equations are partial differential equations with $n$ independent variables, instead of just one independent variable and (2) the multiplier $λ$ inside the differential operator is generally a function of not just the independent variables $\left(x^μ: 0≤μ<n\right)$ but also of the dependent variables $(q,v)$. I don't know if there's even a name for this generalization of the Sturm-Liouville problem.
For such systems, there is a transformation - Classical Renormalization:
$$q → q′ = Zq\quad (Z ≠ 0),$$
with respect to which the form of the Euler-Lagrangian equation is invariant. It can be rewritten as
$$∂·(λ′∂q′) = κ′q′,$$
provided that the equation is rescaled and the coefficients $(κ′,λ′)$ properly defined. Substituting in for $q′$, one gets:
$$∂·(λ′∂(Zq)) = κ′Zq.$$
Rescaling by $Z$, one then has:
$$Z∂·(λ′∂(Zq)) = Z^2κ′q.$$
The left-hand side can be rewritten as:
$$\begin{align}
Z∂·(λ′∂(Zq))
&= \left(Z ∂λ′·∂Z q + Z^2 ∂λ′·∂q\right) + Zλ′ (∂·∂Z q + 2∂Z·∂q + Z ∂·∂q)\\
&= (Z ∂λ′·∂Z + Z λ′∂·∂Z) q + \left(Z^2 ∂λ′ + λ′∂\left(Z^2\right)\right)·∂q + Z^2 λ′∂·∂q\\
&= Z ∂·(λ′∂Z) q + ∂·\left(Z^2λ′∂q\right).
\end{align}$$
Thus, we recover the original equation:
$$∂·\left(Z^2λ′∂q\right) = \left(Z^2κ′ - Z ∂·(λ′∂Z)\right)q,$$
provided we set:
$$κ = Z^2κ′ - Z ∂·(λ′∂Z),\quad λ = Z^2λ′.$$
Inverting these relations, we get:
$$κ′ = \left(\frac{1}{Z}\right)^2 κ - \frac{1}{Z} ∂·\left(λ∂\left(\frac{1}{Z}\right)\right),\quad λ′ = \left(\frac{1}{Z}\right)^2λ,$$
along with
$$q = \frac{1}{Z}q′.$$
The transform is also closed under composition. If we write:
$$q″ = Yq′,\quad κ′ = Y^2κ″ - Y∂·(λ″∂Y),\quad λ′ = Y^2λ″,$$
then it follows that
$$q″ = YZq,\quad κ = (YZ)^2κ″ - YZ∂·(λ″∂(YZ)),\quad λ = (YZ)^2λ″.$$
Therefore, the transformations form a group - the Classical Renormalization Group for the scalar field equations.
The classical theory that is normally quantized has the following properties:
$$κ = -\left(\frac{mc}{ħ}\right)^2λ,\quad ∂_μλ = 0.$$
Seen in this light, the renormalization described above entails not only a rescaling of the coefficients (including the mass $m$), but also a shift of the mass $m$ by $δm$.
In the classical theory that's quantized, the Lagrangian density is "written in by hand" - which may be the very cardinal sin that leads to the problems with infinities! It may be assuming too much to actually write in this or any Lagrangian by hand, and the resulting need to compensate in the back end in the quantized theory with regularization may be the price that's being paid for this.
As yet, I don't have a Rosetta Stone that cross-correlates the Wilsonian approach to what's presented here, but I suspect that this may be what is ultimately underlying the issue.
At the quantum level, a similar exercise can be carried out for the QED Lagrangian density, which can be written (with a gauge-fixing term) as:
$$S = \int 𝔏_Q d^4x,$$
with
$$𝔏_Q = \left(-\frac{ε_0 c}{4} g^{μρ} g^{νσ} F_{μν} F_{ρσ} - \frac{ε_0 c}{2}\left(g^{μν}∂_μA_ν\right)^2 + \bar{ψ} \left(γ^μ (iħ ϝ_2 \overleftrightarrow{∂_μ} + e ϝ_1 A_μ) - ϝ_0\right) ψ\right) \sqrt{|g|},$$
where
$$\overleftrightarrow{∂_μ} ≡ \frac{\overrightarrow{∂_μ} - \overleftarrow{∂_μ}}{2}$$
with the arrows denoting the direction the derivative operator is to be applied in, and where the Lagrangian derivatives for the electromagnetic terms and for the fermion terms are all made explicit and have the relations:
$$ϝ_1 = ϝ_2,\quad ϝ_0 = mc ϝ_1.$$
For electromagnetism, the constitutive coefficient is the reciprocal of the impedance of free space
$$Z = μ_0 c = \frac{1}{ε_0c}$$
up to numeric multiple. Essentially, it's just the permittivity $ε_0$, itself. The field strength is given by
$$F_{μν} = ∂_μA_ν - ∂_νA_μ.$$
Under a renormalization transform, the coefficients become:
$$ϝ_0 = Z_0 ϝ_{0R}, \quad ϝ_1 = Z_1 ϝ_{1R},\quad ϝ_2 = Z_2 ϝ_{2R}, \quad ε_0 = Z_3 ε_R,$$
with
$$Z_1 = Z_2,\quad ϝ_{1R} = ϝ_{2R},\quad ϝ_{0R} = m_Rc ϝ_{1R}.$$
In the usual textbook treatment, the coefficients are suppressed - which basically means that they're mixed in with the field components themselves. So, the respective rescalings are presented as a "bare versus dressed" rescaling on the field components. You can derive what it should be by restoring the original coefficients and moving appropriate combinations of the factors $Z_0$, $Z_1$, $Z_2$ and $Z_3$ into the field components, resulting in the "dressed" fields $ψ_R$, $\bar{ψ}_R$, $A_R$ and $F_R$.
The picture I've presented at the classical level isn't yet complete. If you promote the coefficients to Lagrangian density derivatives, then the corresponding invariants are
$$
I_0 = \bar{ψ}ψ,\quad
I_1 = \bar{ψ} γ^μ A_μ ψ,\quad
I_2 = \bar{ψ} \frac{γ^μ \overrightarrow{∂_μ} - \overleftarrow{∂_μ} γ^μ}2 ψ,\\
I_{3F} = \frac{g^{μρ}g^{νσ}F_{μν}F_{ρσ}}4,\quad
I_{3A} = \frac{\left(g^{μν}∂_μA_ν\right)^2}2.
$$
The corresponding derivatives are:
$$
\frac{∂𝔏_Q}{∂I_0} = -ϝ_0 \sqrt{|g|},\quad
\frac{∂𝔏_Q}{∂I_1} = e ϝ_1 \sqrt{|g|},\quad
\frac{∂𝔏_Q}{∂I_2} = iħ ϝ_2 \sqrt{|g|},\quad
\frac{∂𝔏_Q}{∂I_{3F}} = -ε_0 c \sqrt{|g|} = \frac{∂𝔏_Q}{∂I_{3A}},
$$
and are all assumed to be constant. That may be the "assuming too much" part that leads to the problems that regularization has to fix. (And yes, the Epstein-Glaser / Causal Perturbation Theory people: "distribution splitting" counts as regularization, too!) However, if you promote the coefficients to functions, then the resulting Euler-Lagrange equations - which will (once again) be generalized Sturm-Liouville equations - will not be closed under rescaling! Extra terms are missing that need to be added to get closure, but I'm not sure what or how.
On the quantum side, renormalization group flow is constrained and largely determined by self-consistency; the requirement that counter-terms cancel in the right way. One might be able to determine what the functions for the constitutive coefficients should be, based on this.
