A search for a light pseudoscalar Higgs boson (a) decaying from the 125 GeV (or a heavier) scalar Higgs boson (H) is performed using the 2016 LHC proton-proton collision data at $\sqrt{s} =$ 13 TeV, corresponding to an integrated luminosity of 35.9 fb$^{-1}$, collected by the CMS experiment. The analysis considers gluon fusion and vector boson fusion production of the H, followed by the decay H $\to$ aa $\to$ $\mu\mu\tau\tau$, and considers pseudoscalar masses in the range 3.6 $\lt$ $m_\mathrm{a}$ $\lt$ 21 GeV. Because of the large mass difference between the H and the a bosons and the small masses of the a boson decay products, both the $\mu\mu$ and the $\tau\tau$ pairs have high Lorentz boost and are collimated. The $\tau\tau$ reconstruction efficiency is increased by modifying the standard technique for hadronic $\tau$ lepton decay reconstruction to account for a nearby muon. No significant signal is observed. Model-independent limits are set at 95% confidence level, as a function of $m_\mathrm{a}$, on the branching fraction ($\mathcal{B}$) for H $\to$ aa $\to$ $\mu\mu\tau\tau$, down to 1.5 (2.0) $\times$ 10$^{-4}$ for $m_\mathrm{H} =$ 125 (300) GeV. Model-dependent limits on $\mathcal{B}$(H $\to$ aa) are set within the context of two Higgs doublets plus singlet models, with the most stringent results obtained for Type-III models. These results extend current LHC searches for heavier a bosons that decay to resolved lepton pairs and provide the first such bounds for an H boson with a mass above 125 GeV.