Transient heat transfer plays a major role in many physical and chemical processes, which implies the need for efficient numerical simulation tools in a wide range of engineering fields, such as mechanical, chemical and process engineering. Our aim is to develop a stable, robust, and efficient boundary element algorithm to solve problems arising from applications in those disciplines. It is a well known fact that uniqueness, solvability, and quasi-optimality of space-time Galerkin boundary element methods follows from the boundedness and ellipticity of the thermal single layer- and hyper-singular operator. To handle the arising dense matrix problems we accelerate our solver via the parabolic fast multipole method and thus obtain almost optimal complexity in the number of unknowns. Since the heat kernel is smooth for positive time and exponentially decaying in space, we approximate the kernel by a multivariate Lagrange-Chebyshev interpolation for well separated space-time clusters, which enables an efficient evaluation of their contributions to the thermal layer potentials. Finally, we investigate some benchmark problems to show that the error, induced by the fast multipole method, can be controlled and does not destroy the convergence rates of the Galerkin scheme. Furthermore, we simulate some problems from industrial applications to show the suitability of the proposed method for large scale problems.