  Our aim is to construct high order approximation schemes for general
semigroups of linear operators $P_{t},t\geq 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 )<...<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}}...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,...,r$ such that \[
P_{t}f=\sum_{i=1}^{r}c_{i}\mathbb{E}(P_{t}^{\Pi _{i}(\omega )}f(x))+O(n^{-\nu})
\]% with the expectation concerning the random grids $\Pi _{i}(\omega ).$
Besides, $\text{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.
