Results are presented from a search for the Higgs boson decay , where with or . This search is performed using a sample of proton-proton collision data at a center-of-mass energy of 13~\mathrm{TeV}, recorded by the CMS experiment at the LHC, corresponding to an integrated luminosity of 137~\mathrm{fb}^{-1}. Events are assigned to mutually exclusive categories, which exploit differences in event topology and kinematics of distinct Higgs production modes to enhance signal sensitivity. To detect a potential signal, fits are performed to the distributions of in each of these categories simultaneously. The observed (expected) upper limit at confidence level on the signal strength , defined as the product of the cross section and the branching fraction [\sigma(\mathrm{pp}\to\mathrm{H})\times\,\mathcal{B}(\mathrm{H}\to\mathrm{Z}\gamma)], relative to the standard model expectation, is 4.1\,(). The best fit value of the signal strength is found to be \mu=2.4\,\pm0.9, corresponding to \sigma(\mathrm{pp}\to\mathrm{H})\times\,\mathcal{B}(\mathrm{H}\to\mathrm{Z}\gamma)=0.21\pm0.08~\mathrm{pb}. The statistical significance of the observed excess of events under the background-only hypothesis is 2.7 standard deviations at m_\mathrm{H}=125.38~\mathrm{GeV}.