We show that the running of operators which mix under renormalization can be
computed fully non-perturbatively as a product of continuum step scaling
matrices. These step scaling matrices are obtained by taking the "ratio" of Z
matrices computed at different energies in an RI-MOM type scheme for which
twisted boundary conditions are an essential ingredient. Our method allows us
to relax the bounds of the Rome-Southampton window. We also explain why such a
method is important in view of the light quark physics program of the RBC-UKQCD
collaborations. To illustrate our method, using n_f=2+1 domain-wall fermions,
we compute the non-perturbative running matrix of four-quark operators needed
in K->pipi decay and neutral kaon mixing. Our results are then compared to
perturbation theory.Comment: 8 pages, 7 figures. v2: PRD version, minor changes and few references
adde