This paper proposes a general higher order operator splitting scheme for diffusion semigroups using the Baker-Campbell-Hausdorff type commutator expansion of non-commutative algebra and the Malliavin calculus. An accurate discretization method for the fundamental solution of heat equations or the heat kernel is introduced with a new computational algorithm which will be useful for the inference for diffusion processes. The approximation is regarded as the splitting around the Euler-Maruyama scheme for the density. Numerical examples for diffusion processes are shown to validate the proposed scheme.</jats:p>