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.