You are conflating three conceptually different categories of "regularizations" of seemingly divergent series (and integrals).
The type of resummations that Hardy would talk about are similar to the zeta-function regularization - the example that is most familiar to the physicists. For example,
$$S=\sum_{n=1}^\infty n= -\frac{1}{12}$$
is the most famous sum. Note that this result is unique; it is a well-defined number. In particular, that allows one to calculate the critical dimension of bosonic string theory from $(D-2)S/2+1=0$ and the result is $D=26$. Fundamentally speaking, there is no real divergence in the sum. The "divergent pieces" may be subtracted "completely".
However, in the usual cases of renormalization - of a loop diagram - in quantum field theory, there are divergences. Renormalization removes the "infinite part" of these terms. A finite term is left but the magnitude of the term is not uniquely determined, like it was in the case of the sum of positive integers. Instead, every type of a divergence in a loop diagram produces one parameter - analogous to the coupling constant - that has to be adjusted. Because the finite results can be "anything", this is clearly something else than the zeta-regularization and, more generally, Hardy's procedures whose very goal was to produce unique, well-defined results for seemingly divergent expressions. Infinitesimally speaking, the Renormalization Group only mixes the lower-order contributions (by the number of loops) into a higher-order contribution.
So these are two different things that one should distinguish.
There is another category of problems that is different from both categories above: the summation of the perturbative expansions to all orders. It can be demonstrated that in almost all fields theories - and perturbative string theories as well - the perturbative expansions diverge. For a small coupling, one can sum them up to the smallest term, before the factorial-like coefficient begins to increase the terms again, despite the $g^{2L}$ suppression. The smallest term is of the same order as the leading non-perturbative contributions.
At the very end, if the theory can be non-perturbatively well-defined - and both QCD-like theories and string theory can, at least in principle - the full function as a function of the coupling constant $g$ exists. But it just can't be fully obtained from the perturbative expansion. The Renormalization Group won't really help you because it only mixes the perturbative terms of another order to a perturbative diagram you want to calculate. If you don't know the non-perturbative physics, the equations of the Renormalization Group won't fill the gap because they will keep you in the perturbative realm.
So I have sketched three different things: in the Hardy/zeta problems, the answer to the divergent series was unique; in the particular $L$-loop diagrams in QFT, it wasn't unique but the infinite part was subtracted and the finite part was obtained by a comparison with the experiments; and in the perturbative expansion resummed to all orders, the sum actually didn't converge and indeed, it didn't know about all the information about the full result for a finite $g$.
The last statement may have some subtleties; at least for some theories, the non-perturbative physics is fully determined by the perturbative physics. But I think it is not quite general and we have counterexamples - e.g. for AdS/CFT with orthogonal groups and different discrete values of $B$ etc. So it means that the perturbative expansion doesn't uniquely determine the theory non-perturbatively.
Because the three examples differ at the level of "what can be calculated" and "what cannot", they are different.
At least at the operational level, whether an operator is relevant or irrelevant (in the IR) tells you about it's canonical scaling dimension.
I think the BPHZ picture of renormalization might help here. Given a physical theory, we'd like to estimate which Feynman graphs of that theory diverge. You can estimate a "superficial degree of divergence" for each Feynman diagram by considering the canonical scaling dimension of each vertex or edge in the diagram. Say you have a diagram with a vertex corresponding to an IR-irrelevant (UV-relevant) interaction term. Then, in the UV (for large momenta), it will cause the amplitude to diverge. Now, as you go to higher loop order, or insert more of those vertices at same loop order, the divergence becomes worse. And as you add more and more such vertices, the corresponding diagrams become more and more important in the UV, due to the couplings being UV-relevant! The essence is that there are infinitely many ("independent") Feynman diagrams which diverge, so you can't cook up enough counterterms to cancel all these divergences. So these terms in the lagrangian give rise to non-renormalizable theories.
See, for eg. Peskin & Schroeder section 10.1, or the link above for how to calculate the superficial degree of divergence of any diagram. Another technicality is that you don't actually consider Feynman graphs, you instead consider subgraphs. What that means is that the external (uncontracted) legs can be off-shell, as if they were sitting inside a bigger Feynman diagram.
Given your question, maybe you've already come across the stuff I've said and are asking for something else, in which case your question is not clear to me.
Best Answer
There are several interesting questions in the main question, plus a point in the comment that I want to address. Similar ideas are discussed here and in arxiv 0702.365.
Disclaimer : I will only speak about QFT that have a finite UV cut-off $\Lambda$. That remove all the complications of the definition of the (non-pertubative) continuum limit $\Lambda\to \infty$, as discussed in the answer linked above. It is only if you believe that a particular QFT is the absolute, true description of the universe that this limit is interesting. And it's pretty sure that it is not the case. Nevertheless, we can be interested in the limit where $\Lambda$ is very large compare to all other energy scales (which is equivalent to take $\Lambda\to\infty$).
First of all, bare parameters can be physical, and measured. They just don't correspond to the same quantities than renormalized parameters. For instance, take the (classical) Ising model. It has one coupling constant $K=J/T$. Using standard calculations, one can rewrite the partition function as a field theory with action $$ S[\phi_i]=\sum_{ij} t_{ij} \phi_i \phi_j- \sum_i \ln \cosh \phi_i, $$ where $\phi_i$ is the value of the field on a lattice site $i$, and $t_{ij}$ is related to the interaction energy $K$ (see for example this article for the details). Thus, if you know $K$, which is accessible (for instance in the case of simulations), you know the bare parameters ! You can also measure this kind of parameters experimentally (exchange energy).
Notice that this field theory is non-perturbative (if one expands the potential, all the coupling constants (and there are an infinity of them) are of the same order ! The only reason one can use the pertubative $\phi^4$ theory to describe the Ising model is because one is usually interested only in the universal quantities, that don't care about the details of the microscopic theory, as long as the universality class is the same.
For the sake of completeness : if one expand the potential (the $\ln\cosh$) to the order four, the quadratic term will be called the mass term and the fourth term the interaction. There is also a contribution to the mass coming from $t_{ii}$. One can then show that for $K$ large enough, the potential has two non trivial minimum, corresponding to the ferromagnetic phase. The critical value of $K$ for the transition, noted $K^0_c$, is at the mean-field level both wrong and non-physical.
So far so good. Now let's add the fluctuations of the field, that will "renormalize" the theory. At one loop, one sees that the quadratic term (the "mass") has a correction which is proportional to some power of the cut-off, that is, it depends on the way the regulation is made, whether the lattice is cubic or triangular, etc. Is it a problem ? Not at all. This is just telling you that the (real, physical) critical coupling $K_c$ is non-universal, it depends strongly on the microscopic details of the system. In some sense, the calculation of these "divergent" (in reality, cut-off dependent) integrals corresponds to the calculation of the critical temperature, knowing the microscopic physics.
From the Wilsonian RG, we know that some quantities will depend on the cut-off, such as the critical temperature. They are usually very hard to compute using a field theory, since the initial action is non-pertubative and one can not use the pertubative RG (wilsonian or not). Only non-pertubative schemes (numerical approaches, or the non-pertubative RG discussed in the arxiv articles linked above) can access these quantities. But there are universal quantities, such as the critical exponents, that can be computed with pertubative approach, as long as one stays in the same universality class. To compute these quantities, one has to be close to the fixed point of the RG, that is, at energy very low compare to the microscopic scales, equivalent to take $\Lambda\to \infty$.
This leads us to see the difference between Wilson RG and "old school" RG. In the former, one is imposing the microscopic value of $K$, and then look at what is the physical mass. In the latter, one imposes the physical value of the mass, and does not care about the microscopic details, so one wants to send $\Lambda\to \infty$. One thus has to absorb the "correction" to the mass in order to fix it.
So, to (finally, but partially) answer your question "The running of the coupling in Wilson's approach has nothing to do with the bare parameters going to infinity when the cutoff is removed right?" :
In the Wilsonian approach, one starts from the microscopic scale $\Lambda$ and looks at what's going at smaller energy, whereas in the "standard" approach, one fixes than macroscopic scale and sends $\Lambda\to \infty$ in order to effectively probe smaller and smaller energy scales.