  We present an improved result of lattice computation of the proton decay
matrix elements in $N_f=2+1$ QCD. In this study, the significant improvement of
statistical accuracy by adopting the error reduction technique of
All-mode-averaging, is achieved for relevant form factor to proton (and also
neutron) decay on the gauge ensemble of $N_f=2+1$ domain-wall fermions in
$m_\pi=0.34$--0.69 GeV on 2.7~fm$^3$ lattice as used in our previous work
\cite{Aoki:2013yxa}. We improve total accuracy of matrix elements to 10--15\%
from 30--40\% for $p\rightarrow\pi e^+$ or from 20--40\% for $p\rightarrow K
\bar\nu$. The accuracy of the low energy constants $\alpha$ and $\beta$ in the
leading-order baryon chiral perturbation theory (BChPT) of proton decay are
also improved. The relevant form factors of $p\rightarrow \pi$ estimated
through the "direct" lattice calculation from three-point function appear to be
1.4 times smaller than those from the "indirect" method using BChPT with
$\alpha$ and $\beta$. It turns out that the utilization of our result will
provide a factor 2--3 larger proton partial lifetime than that obtained using
BChPT. We also discuss the use of these parameters in a dark matter model.
