Our aim is to construct high order approximation schemes for general semigroups of linear operators \(P_{t},t\ge 0\). In order to do it, we fix a time horizon T and the discretization steps \(h_{l}=\frac{T}{n^{l}},l\in {\mathbb {N}}\) and we suppose that we have at hand some short time approximation operators \(Q_{l}\) such that \(P_{h_{l}}=Q_{l}+O(h_{l}^{1+\alpha })\) for some \(\alpha >0\). Then, we consider random time grids \(\Pi (\omega )=\{t_0(\omega )=0<t_{1}(\omega )<\cdots <t_{m}(\omega )=T\}\) such that for all \(1\le k\le m\), \(t_{k}(\omega )-t_{k-1}(\omega )=h_{l_{k}}\) for some \(l_{k}\in {\mathbb {N}}\), and we associate the approximation discrete semigroup \(P_{T}^{\Pi (\omega )}=Q_{l_{n}}\ldots Q_{l_{1}}\). Our main result is the following: for any approximation order \(\nu \), we can construct random grids \(\Pi _{i}(\omega )\) and coefficients \(c_{i}\), with \(i=1,\ldots ,r\) such that $$\begin{aligned} P_{t}f=\sum _{i=1}^{r}c_{i}{\mathbb {E}}(P_{t}^{\Pi _{i}(\omega )}f(x))+O(n^{-\nu }) \end{aligned}$$with the expectation concerning the random grids \(\Pi _{i}(\omega ).\) Besides, \(Card (\Pi _{i}(\omega ))=O(n)\) and the complexity of the algorithm is of order n, for any order of approximation \(\nu \). The standard example concerns diffusion processes, using the Euler approximation for \(Q_l\). In this particular case and under suitable conditions, we are able to gather the terms in order to produce an estimator of \(P_tf\) with finite variance. However, an important feature of our approach is its universality in the sense that it works for every general semigroup \(P_{t}\) and approximations. Besides, approximation schemes sharing the same \(\alpha \) lead to the same random grids \(\Pi _{i}\) and coefficients \(c_{i}\). Numerical illustrations are given for ordinary differential equations, piecewise deterministic Markov processes and diffusions.