We investigate the prospects for the discovery of neutral Higgs bosons with a pair of muons by direct searches at the CERN Large Hadron Collider (LHC) as well as by indirect searches in the rare decay $B_s to mu^+mu^-$ at the Fermilab Tevatron and the LHC. Promising results are found for the minimal supersymmetric standard model, the minimal supergravity (mSUGRA) model, and supergravity models with non-universal Higgs masses (NUHM SUGRA). For $tanbeta simeq 50$, we find that (i) the contours for a branching fraction of $B(B_s to mu^+mu^-) = 1 times 10^{-8}$ in the parameter space are very close to the $5sigma$ contours for $pp to bphi^0 to bmu^+mu^- +X, phi^0 = h^0, H^0, A^0$ at the LHC with an integrated luminosity ($L$) of 30 fb$^{-1}$, (ii) the regions covered by $B(B_s to mu^+mu^-) ge 5times 10^{-9}$ and the discovery region for $bphi^0 to bmu^+mu^-$ with 300 fb$^{-1}$ are complementary in the mSUGRA parameter space, (iii) in NUHM SUGRA models, a discovery of $B(B_s to mu^+mu^-) simeq 5times 10^{-9}$ at the LHC will cover regions of the parameter space beyond the direct search for $bphi^0 to bmu^+mu^-$ with $L = 300$ fb$^{-1}$.