  The Picard-Lefschetz theory has been attracting much attention as a tool to
evaluate a multi-variable integral with a complex weight, which appears in
various important problems in theoretical physics. The idea is to deform the
integration contour based on Cauchy's theorem using the so-called gradient flow
equation. In this paper, we propose a fast Hybrid Monte Carlo algorithm for
evaluating the integral, where we "backpropagate" the force of the fictitious
Hamilton dynamics on the deformed contour to that on the original contour,
thereby reducing the required computational cost by a factor of the system
size. Our algorithm can be readily extended to the case in which one integrates
over the flow time in order to solve not only the sign problem but also the
ergodicity problem that occurs when there are more than one thimbles
contributing to the integral. This enables, in particular, efficient
identification of all the dominant saddle points and the associated thimbles.
We test our algorithm by calculating the real-time evolution of the wave
function using the path integral formalism.
