A search is presented for a heavy resonance $Y$ decaying into a Standard Model Higgs boson $H$ and a new particle $X$ in a fully hadronic final state. The full Large Hadron Collider Run 2 dataset of proton-proton collisions at $\sqrt{s}= 13$ TeV collected by the ATLAS detector from 2015 to 2018 is used, and corresponds to an integrated luminosity of 139 fb$^{-1}$. The search targets the high $Y$-mass region, where the $H$ and $X$ have a significant Lorentz boost in the laboratory frame. A novel signal region is implemented using anomaly detection, where events are selected solely because of their incompatibility with a learned background-only model. It is defined using a jet-level tagger for signal-model-independent selection of the boosted $X$ particle, representing the first application of fully unsupervised machine learning to an ATLAS analysis. Two additional signal regions are implemented to target a benchmark $X$ decay into two quarks, covering topologies where the $X$ is reconstructed as either a single large-radius jet or two small-radius jets. The analysis selects Higgs boson decays into $b\bar{b}$, and a dedicated neural-network-based tagger provides sensitivity to the boosted heavy-flavor topology. No significant excess of data over the expected background is observed, and the results are presented as upper limits on the production cross section $\sigma(pp \rightarrow Y \rightarrow XH \rightarrow q\bar{q}b\bar{b}$) for signals with $m_Y$ between 1.5 and 6 TeV and $m_X$ between 65 and 3000 GeV.