A search for the non-resonant production of Higgs boson pairs (HH) via the gluon-gluon and vector boson fusion processes in final states with two bottom quarks and two tau leptons is presented. The search uses data from proton-proton collisions at a center-of-mass energy of $\sqrt{s}=13~\mathrm{TeV}$ recorded with the CMS detector at the LHC, corresponding to an integrated luminosity of $138~\mathrm{fb}^{-1}$. Events in which at least one tau lepton decays hadronically are considered and multiple machine learning techniques are used to identify and extract the signal. The data are found to be consistent, within uncertainties, with the standard model (SM) background predictions. Upper limits on the HH production cross section and constraints on anomalous Higgs boson couplings are set. The observed (expected) upper limit at 95% confidence level corresponds to 3.3 (5.2) times the SM prediction for the inclusive HH cross section and to 124 (154) times the SM prediction for the vector boson fusion HH cross section. At a 95% confidence level, the Higgs field self-coupling is constrained to be within -1.8 and 8.8 times the standard model expectation, and the coupling of two Higgs bosons to two vector bosons is constrained to be within -0.4 and 2.6.