A search for new heavy resonances decaying to WW, WZ, ZZ, WH, or ZH boson pairs in the all-jets final state is presented. The analysis is based on proton-proton collision data recorded by the CMS detector in 2016-2018 at a centre-of-mass energy of 13 TeV at the CERN LHC, corresponding to an integrated luminosity of 138 fb$^{-1}$. The search is sensitive to resonances with masses between 1.3 and 6 TeV, decaying to bosons that are highly Lorentz-boosted such that each of the bosons forms a single large-radius jet. Machine learning techniques are employed to identify such jets. No significant excess over the estimated standard model background is observed. A maximum local significance of 3.6 standard deviations, corresponding to a global significance of 2.3 standard deviations, is observed at masses of 2.1 and 2.9 TeV. In a heavy vector triplet model, spin-1 Z' and W' resonances with masses below 4.8 TeV are excluded at the 95% confidence level (CL). These limits are the most stringent to date. In a bulk graviton model, spin-2 gravitons and spin-0 radions with masses below 1.4 and 2.7 TeV, respectively, are excluded at 95% CL. Production of heavy resonances through vector boson fusion is constrained with upper cross section limits at 95% CL as low as 0.1 fb.