This work presents a new physically-motivated supervised machine learning method, Hydro-BAM, to reproduce the three-dimensional Lyman-$alpha$ forest field in real and in redshift space learning from a reference hydrodynamic simulation, thereby saving about 7 orders of magnitude in computing time. We show that our method is accurate up to $ksim1,h,rm{Mpc}^{-1}$ in the one- (PDF), two- (power-spectra) and three-point (bi-spectra) statistics of the reconstructed fields. When compared to the reference simulation including redshift space distortions, our method achieves deviations of $lesssim2%$ up to $k=0.6,h,rm{Mpc}^{-1}$ in the monopole, $lesssim5%$ up to $k=0.9,h,rm{Mpc}^{-1}$ in the quadrupole. The bi-spectrum is well reproduced for triangle configurations with sides up to $k=0.8,h,rm{Mpc}^{-1}$. In contrast, the commonly-adopted Fluctuating Gunn-Peterson approximation shows significant deviations already neglecting peculiar motions at configurations with sides of $k=0.2-0.4,h,rm{Mpc}^{-1}$ in the bi-spectrum, being also significantly less accurate in the power-spectrum (within 5$%$ up to $k=0.7,h,rm{Mpc}^{-1}$). We conclude that an accurate analysis of the Lyman-$alpha$ forest requires considering the complex baryonic thermodynamical large-scale structure relations. Our hierarchical domain specific machine learning method can efficiently exploit this and is ready to generate accurate Lyman-$alpha$ forest mock catalogues covering large volumes required by surveys such as DESI and WEAVE.