Summary: Projection methods have emerged as competitive techniques for solving large scale matrix Lyapunov equations. We explore the numerical solution of this class of linear matrix equations when a minimal residual (MR) condition is used during the projection step. We derive both a new direct method, and a preconditioned operator-oriented iterative solver based on a conjugate gradient least squares method, for solving the projected reduced least squares problem. Numerical experiments with benchmark problems show the effectiveness of an MR approach over a Galerkin procedure using the same approximation space.