The limited-memory quasi-Newton method with simple bounds is used to develop a novel, fully 3D magnetotelluric (MT) inversion technique. This nonlinear inversion is based on iterative minimization of a classical Tikhonov regularized penalty function. However, instead of the usual model space of log resistivities, the approach iterates in a model space with simple bounds imposed on the conductivities of the 3D target. The method requires storage proportional to , where N is the number of conductivities to be recovered and is the number of correction pairs (practically, only a few). These requirements are much less than those imposed by other Newton methods, which usually require storage proportional to or , where M is the number of data to be inverted. The derivatives of the penalty function are calculated using an adjoint method based on electromagnetic field reciprocity. The inversion involves all four entries of the MT impedance matrix; the integral equation forward-modeling code is used as an engine for this inversion. Convergence, performance, and accuracy of the inversion are demonstrated on synthetic numerical examples. After investigating erratic resistivities in the upper part of the model obtained for one of the examples, we conclude that the standard Tikhonov regularization is not enough to provide consistently smooth underground structures. An additional regularization helps to overcome the problem.