  We formulate the path integral of two- and three-flavor Wilson fermion in two
dimensions as a multilayer Grassmann tensor network by the matrix product
decomposition. Thanks to this new description, the memory cost scaling is
reduced from $\mathrm{O}(\mathrm{e}^{N_{f}})$ for the conventional construction
to $\mathrm{O}(N_{f})$. Based on this representation, we develop a
coarse-graining algorithm where spatially or temporally adjacent Grassmann
tensors are converted into a canonical form along a virtual direction before we
carry out the spacetime coarse-graining. Benchmarking with the lattice
Gross-Neveu model at finite density, we see that the Silver Blaze phenomenon in
the pressure and number density is captured with relatively small bond
dimensions.
