A recently developed variant of the so-called optimized perturbation theory (OPT), making it perturbatively consistent with renormalization group (RG) properties, RGOPT, was shown to drastically improve its convergence for zero temperature theories. Here the RGOPT adapted to finite temperature is illustrated with a detailed evaluation of the two-loop pressure for the thermal scalar $ \lambda\phi^4$ field theory. We show that already at the simple one-loop level this quantity is exactly scale-invariant by construction and turns out to qualitatively reproduce, with a rather simple procedure, results from more sophisticated resummation methods at two-loop order, such as the two-particle irreducible approach typically. This lowest order also reproduces the exact large-$N$ results of the $O(N)$ model. Although very close in spirit, our RGOPT method and corresponding results differ drastically from similar variational approaches, such as the screened perturbation theory or its QCD-version, the (resummed) hard thermal loop perturbation theory. The latter approaches exhibit a sensibly degrading scale dependence at higher orders, which we identify as a consequence of missing RG invariance. In contrast RGOPT gives a considerably reduced scale dependence at two-loop level, even for relatively large coupling values $\sqrt{\lambda/24}\sim {\cal O}(1)$, making results much more stable as compared with standard perturbation theory, with expected similar properties for thermal QCD.