  We apply non-perturbative renormalization to bilinears composed of improved
staggered fermions. We explain how to generalize the method to staggered
fermions in a way which is consistent with the lattice symmetries, and
introduce a new type of lattice bilinear which transforms covariantly and
avoids mixing. We derive the consequences of lattice symmetries for the
propagator and vertices. We implement the method numerically for
hypercubic-smeared (HYP) and asqtad valence fermion actions, using lattices
with asqtad sea quarks generated by the MILC collaboration. We compare the
non-perturbative results so obtained to those from perturbation theory, using
both scale-independent ratios of bilinears (of which we calculate 26), and the
scale-dependent bilinears themselves. Overall, we find that one-loop
perturbation theory provides a successful description of the results for
HYP-fermions if we allow for a truncation error of roughly the size of the
square of the one-loop term (for ratios) or of size O(1) \times \alpha^2 (for
the bilinears themselves). Perturbation theory is, however, less successful at
describing the non-perturbative asqtad results.
