  Hadronic matrix elements of proton decay are essential ingredients to bridge
the grand unification theory to low energy observables like proton lifetime. In
this paper we non-perturbatively calculate the matrix elements, relevant for
the process of a nucleon decaying into a pseudoscalar meson and an anti-lepton
through generic baryon number violating four-fermi operators. Lattice QCD with
2+1 flavor dynamical domain-wall fermions with the {\it direct} method, which
is direct measurement of matrix element from three-point function without
chiral perturbation theory, are used for this study to have good control over
the lattice discretization error, operator renormalization, and chiral
extrapolation. The relevant form factors for possible transition process from
an initial proton or neutron to a final pion or kaon induced by all types of
three quark operators are obtained through three-point functions of
(nucleon)-(three-quark operator)-(meson) with physical kinematics. In this
study all the relevant systematic uncertainties of the form factors are taken
into account for the first time, and the total error is found to be the range
30%-40% for $\pi$ and 20%-40% for $K$ final states.
