  We explore double Higgs production via gluon fusion in the $b\bar{b} \gamma
\gamma $ channel at the high-luminosity LHC using machine learning tools. We
first propose a Bayesian optimization approach to select cuts on kinematic
variables, obtaining a $30-50$ \% increase in the significance compared to
current results in the literature. We show that this improvement persists once
systematic uncertainties are taken into account. We next use boosted decision
trees (BDT) to further discriminate signal and background events. Our analysis
shows that a joint optimization of kinematic cuts and BDT hyperparameters
results in an appreciable improvement in the significance. Finally, we perform
a multivariate analysis of the output scores of the BDT. We find that assuming
a very low level of systematics, the techniques proposed here will be able to
confirm the production of a pair of Standard Model Higgs bosons at 5$\sigma$
level with 3 ab$^{-1}$ of data. Assuming a more realistic projection of the
level of systematics, around 10\%, the optimization of cuts to train BDTs
combined with a multivariate analysis delivers a respectable significance of
4.6$\sigma$. Even assuming large systematics of 20\%, our analysis predicts a
3.6$\sigma$ significance, which represents at least strong evidence in favor of
double Higgs production. We carefully incorporate background contributions
coming from light flavor jets or $c$-jets being misidentified as $b$-jets and
jets being misidentified as photons in our analysis.
