A search for new heavy resonances decaying to a pair of Higgs bosons in proton-proton collisions at a center-of-mass energy of 13 TeV is presented. Data were collected by the CMS detector at the LHC from 2016 to 2018, corresponding to an integrated luminosity of 138 . The search considers resonances with a mass between 0.8 and 4.5 TeV using events in which one Higgs boson decays into a bottom quark-antiquark pair and the other decays into final states with either one or two charged leptons. Specifically, these include the single-lepton final state of the decay and the dilepton final states of both the and decays, where in the final state corresponds to or . The signal is extracted using a two-dimensional maximum likelihood fit of the jet mass and invariant mass distributions. No significant excess above the standard model expectation is observed in data. Model-independent exclusion limits are placed on the the product of the cross section and branching fraction () for spin-0 and spin-2 massive bosons decaying to . The results are interpreted in the context of radion and bulk graviton production in models with a warped extra spatial dimension.