  For the first time, we present a lattice QCD determination of Mellin moments
of unpolarized generalized parton distributions (GPDs) of the proton from an
analysis of the quasi-GPD matrix elements within the short-distance
factorization framework. We perform our calculation on an $N_f$=2+1+1 twisted
mass fermions ensemble with a clover improvement at lattice spacing $a=0.093$
fm and a pion mass of $m_\pi=260$ MeV. Focusing on the zero-skewness case, the
iso-vector and iso-scalar quasi-GPDs are calculated from the $\gamma_0$
definition, as well as a recently proposed Lorentz-invariant definition. We
utilize data on both symmetric and asymmetric kinematic frames, which allows us
to obtain the Mellin moments for several values of the momentum transfer, $-t$,
in the range 0.17 to $2.77~\rm{GeV}^2$. We use the ratio scheme for GPDs, i.e.
renormalization group invariant ratios with leading-twist factorization formula
and perturbatively calculated matching coefficients up to the
next-next-to-leading order (NNLO) to extract Mellin moments of GPDs, which are
consistent with renormalization-group improved results. We compare our
determination from quasi-GPDs with the results extracted using standard
calculations of Mellin moments of local operators, specifically those related
to the electromagnetic and gravitational form factors. We estimated the moments
of GPDs up to the fifth ones for the first time. By extrapolating the Mellin
moments to $-t=0$, we obtained the quark charges, momentum fraction, as well as
the angular momentum contributions to the proton spin. The impact parameter
space interpretation of the GPD moments is discussed, which provides insights
into the spatial distribution of unpolarized quarks and their correlations in
the transverse plane of an unpolarized or transversely polarized proton.
