We consider the numerical solution of linear partial differential equations with coefficients that vary on scales much smaller than the computational domain. Naturally, a direct discretization on a grid coarser than those scales will in general give an inaccurate discrete approximation. Instead, one can derive an effective coarse grid operator whose structure is similar to the one given by direct discretization, but with a locally altered stencil that takes the effect of subgrid scales into account. In analogy with classical homogenization, we call this operator the numerically homogenized operator. We discuss a general procedure for deriving such operators. It is based on wavelet projections of the discrete fine grid operator followed by sparse approximation. We review some theoretical results related to this method and show a few numerical examples.