We present a new determination of the $B_s$ leptonic decay constant from lattice QCD simulations that use gluon configurations from MILC and a highly improved discretization of the relativistic quark action for both valence quarks. Our result, $f_{B_s} = 0.225(4)$,GeV, is almost three times more accurate than previous determinations. We analyze the dependence of the decay constant on the heavy quarks mass and obtain the first empirical evidence for the leading $1/sqrt{m_h}$ dependence predicted by Heavy Quark Effective Theory (HQET). As a check, we use our analysis technique to calculate the $m_{B_s}-m_{eta_b}/2$ mass difference. Our result agrees with experiment to within errors of $11,mathrm{MeV}$ (better than 2%). We discuss how to extend our analysis to other quantities in $B_s$ and $B$ physics, making 2%-precision possible for the first time.