In this paper, we investigate the HAL QCD potential in the $I=1$ $pi pi$ scattering using the hybrid method for all-to-all propagators, in which a propagator is approximated by low-eigenmodes and the remaining high-eigenmode part is stochastically estimated. To verify the applicability of the hybrid method to systems containing quark creation$/$annihilation contributions such as the $rho$ meson, we calculate the $I=1$ $pipi$ potential with the 2+1 flavor gauge configurations on $16^3 times 32$ lattice at the lattice spacing $a approx 0.12$ fm and $(m_{pi},m_{rho}) approx (870, 1230)$ MeV, in which the $rho$ meson appears as a deeply-bound state. While we find that the naive stochastic evaluations for quark creation$/$annihilation contributions lead to extremely large statistical fluctuations, additional noise reduction methods enable us to obtain a sufficiently precise potential, which shows a strong attractive force. We also confirm that the binding energy and $k^3 cot delta$ obtained from our potential are roughly consistent with an existing $rho$ meson bound state, within a large systematic error associated with our calculation, whose possible origin is also discussed.