A long-standing issue in peculiar velocity cosmology is whether the halo/galaxy velocity bias b _ { v } = 1 at large scale . The resolution of this important issue must resort to high precision cosmological simulations . However , this is hampered by another long-standing “ sampling artifact ” problem in volume weighted velocity measurement . We circumvent this problem with a hybrid approach . We first measure statistics free of sampling artifact , then link them to volume weighted statistics in theory , finally solve for the velocity bias . b _ { v } determined by our method is not only free of sampling artifact , but also free of cosmic variance . We apply this method to a \Lambda CDM N-body simulation of 3072 ^ { 3 } particles and 1200 { Mpc } / h box size . For the first time , we determine the halo velocity bias to 0.1 \% - 1 \% accuracy . Our major findings are as follows : ( 1 ) b _ { v } \neq 1 at k > 0.1 h / { Mpc } . The deviation from unity ( |b _ { v } -1 | ) increases with k . Depending on halo mass and redshift , it may reach \mathcal { O } ( 0.01 ) at k = 0.2 h / { Mpc } and \mathcal { O } ( 0.05 ) at k \sim 0.3 h / { Mpc } . The discovered b _ { v } \neq 1 has statistically significant impact on structure growth rate measurement by spectroscopic redshift surveys , including DESI , Euclid and SKA . ( 2 ) Both the sign and the amplitude of b _ { v } -1 depend on mass and redshift . These results disagree with the peak model prediction in that b _ { v } has much weaker deviation from unity , varies with redshift , and can be bigger than unity . ( 3 ) Most of the mass and redshift dependences can be compressed into a single dependence on the halo density bias . Based on this finding , we provide an approximate two-parameter fitting formula .