A search for pairs of Higgs bosons produced via gluon and vector boson fusion is presented, focusing on the four $b$ quark final state. The data sample consists of proton-proton collisions at a center-of-mass energy of 13 TeV, collected with the CMS detector at the LHC, and corresponds to an integrated luminosity of $138\text{ }\text{ }{\mathrm{fb}}^{\ensuremath{-}1}$. No deviation from the background-only hypothesis is observed. A 95% confidence level upper limit on the Higgs boson pair production cross section is observed at 3.9 times the standard model prediction for an expected value of 7.8. Constraints are also set on the modifiers of the Higgs field self-coupling, ${\ensuremath{\kappa}}_{\ensuremath{\lambda}}$, and of the coupling of two Higgs bosons to two vector bosons, ${\ensuremath{\kappa}}_{2\text{ }\mathrm{V}}$. The observed (expected) allowed intervals at the 95% confidence level are $\ensuremath{-}2.3<{\ensuremath{\kappa}}_{\ensuremath{\lambda}}<9.4$ ($\ensuremath{-}5.0<{\ensuremath{\kappa}}_{\ensuremath{\lambda}}<12.0$) and $\ensuremath{-}0.1<{\ensuremath{\kappa}}_{2\text{ }\mathrm{V}}<2.2$ ($\ensuremath{-}0.4<{\ensuremath{\kappa}}_{2\text{ }\mathrm{V}}<2.5$). These are the most stringent observed constraints to date on the $HH$ production cross section and on the ${\ensuremath{\kappa}}_{2\text{ }\mathrm{V}\text{ }}$ coupling.