  We calculate non-perturbative renormalization factors at hadronic scale for
$\Delta S=2$ four-quark operators in quenched domain-wall QCD using the
Schr\"{o}dinger functional method. Combining them with the non-perturbative
renormalization group running by the Alpha collaboration, our result yields the
fully non-perturbative renormalization factor, which converts the lattice bare
$B_K$ to the renormalization group invariant (RGI) $\hat{B}_K$. Applying this
to the bare $B_K$ previously obtained by the CP-PACS collaboration at
$a^{-1}\simeq 2, 3, 4$ GeV, we obtain $\hat{B}_K=0.782(5)(7)$ (equivalent to
$B_K^{\bar{\rm MS}}({\rm NDR}, 2 {\rm GeV}) = 0.565(4)(5)$ by 2-loop running)
in the continuum limit, where the first error is statistical and the second is
systematic due to the continuum extrapolation. Except the quenching error, the
total error we have achieved is less than 2%, which is much smaller than the
previous ones. Taking the same procedure, we obtain $m_{u,d}^{\rm
RGI}=5.613(66)$ MeV and $m_s^{\rm RGI}=147.1(17)$ MeV (equivalent to
$m_{u,d}^{\bar{\rm MS}}(2 {\rm GeV})=4.026(48)$ MeV and $m_{s}^{\bar{\rm MS}}(2
{\rm GeV})=105.6(12)$ MeV by 4-loop running) in the continuum limit.
