Blue horizontal branch stars (BHBs), excellent distant tracers for probing the Milky Way’s halo density profile, are distinguished in the g−r0 versus (i − z)0 color space from another class of stars, blue straggler stars. We develop a Bayesian mixture model to classify BHBs using high-precision photometry data from the Dark Energy Survey Data Release 2 (DES DR2). We select ∼2100 highly probable BHBs based on their griz photometry and the associated uncertainties, and we use these stars to map the stellar halo over the Galactocentric radial range 20 kpc ≲ R ≲ 70 kpc. After excluding known stellar overdensities, we find that the number density n ⋆ of BHBs can be represented by a power-law density profile n ⋆ ∝ R −α with an index of α=4.34−0.12+0.13±0.52 , consistent with existing literature values. In addition, we examine the impact of systematic errors and the spatial inhomogeneity on the fitted density profile. Our work demonstrates the effectiveness of high-precision griz photometry in selecting BHBs. The upcoming photometric survey from the Rubin Observatory, expected to reach depths 2–3 mag greater than DES during its 10 yr mission, will enable us to investigate the density profile of the Milky Way’s halo out to the virial radius, unraveling the complex processes of formation and evolution in our Galaxy.