We develop a novel approach in exploring the joint dependence of halo bias on multiple halo properties using Gaussian process regression .
Using a \Lambda CDM N -body simulation , we carry out a comprehensive study of the joint bias dependence on halo structure , formation history and environment .
We show that the bias is a multivariate function of halo properties that falls into three regimes .
For massive haloes , halo mass explains the majority of bias variation .
For early-forming haloes , bias depends sensitively on the recent mass accretion history .
For low-mass and late-forming haloes , bias depends more on the structure of a halo such as its shape and spin .
Our framework enables us to convincingly prove that V _ { max } / V _ { vir } is a lossy proxy of formation time for bias modelling , whereas the mass , spin , shape and formation time variables are non-redundant with respect to each other .
Combining mass and formation time largely accounts for the mass accretion history dependence of bias .
Combining all the internal halo properties fully accounts for the density profile dependence inside haloes , and predicts the clustering variation of individual haloes to a 20 \% level at \sim 10 { Mpc } h ^ { -1 } .
When an environmental density is measured outside 1 { Mpc } h ^ { -1 } from the halo centre , it outperforms and largely accounts for the bias dependence on the internal halo structure , explaining the bias variation above a level of 30 \% .