The solar atmosphere is diagnosed by solving the polarized radiative transfer problem for plasmas in Non-Local Thermodynamical Equilibrium (NLTE). A key challenge in multidimensional NLTE diagnosis is to integrate the radiative transfer equation (RTE) efficiently, but current methods are local, limited to constant propagation matrices. The formalism presented in this paper lays the foundations to achieve an efficient non-local integration of the RTE based on the Magnus expansion. Our approach starts formulating the problem in terms of rotations represented in the Lorentz-Poincare group (Stokes formalism) and motivating the use of the Magnus expansion. Then, we combine a highly detailed algebraic characterization of the propagation matrix with Magnus to reformulate the homogenous solution to the RTE. Thus, we obtain a compact analytical evolution operator supporting arbitrary variations of the propagation matrix to first and second order in the Magnus expansion and showing the way to higher orders. Finally, we reformulate the inhomogeneous part to make it solvable with the Magnus expansion too, which gives rise to an interesting new object: the inhomogeneous evolution operator. This provides the first consistent and general formal solution of the RTE that furthermore is non-local, efficient, and accurate by design, with other peculiarity: it separates the integration part from the algebraic formal solution. Our framework is verified analytically and with a computational implementation that leads to a new family of numerical radiative transfer methods and suggests several applications, e.g. accelerating NLTE calculations with non-local radiative transfer. With cosmetic changes, our results apply to other universal physical problems sharing the Lorentz-Poincare algebra of the RTE and special relativity.