  The recently developed stochastic gradient method combined with Monte Carlo
sampling techniques [PRB {\bf 95}, 195154 (2017)] offers a low scaling and
accurate method to optimize the projected entangled pair states (PEPS). We
extended this method to the fermionic PEPS (fPEPS). To simplify the
implementation, we introduce a fermi arrow notation to specify the order of the
fermion operators in the virtual entangled EPR pairs. By defining some local
operation rules associated with the fermi arrows, one can implement fPEPS
algorithms very similar to that of standard PEPS. We benchmark the method for
the interacting spinless fermion models, and the t-J models. The numerical
calculations show that the gradient optimization greatly improves the results
of simple update method. Furthermore, much larger virtual bond dimensions ($D$)
and truncation dimensions ($D_c$) than those of boson and spin systems are
necessary to converge the results. The method therefore offer a powerful tool
to simulate fermion systems because it has much lower scaling than the direct
contraction methods.
