We present an inversion method for 3D electrical imaging in media with an inhomogeneous and anisotropic conductivity distribution. The conductivity distribution is discretized via finite elements and is described by a second-order tensor at each finite element node. The inversion method is formulated as a functional optimization with an error functional containing terms measuring data misfit and model covariance by means of smoothness, anisotropy and deviation from a starting model. Including the model covariance information overcomes the problem of ill-posedness at the expense of limiting the allowed models to the class of models which are compatible with the provided model covariance information. The discretized form of the error functional is minimized by a Levenberg–Marquardt type method using an iterative preconditioned conjugate gradient solver. The use of an iterative solver allows one to bypass the actual computation of the Jacobian or an inverse system matrix. The use of a memory efficient iterative solver together with the implementation on parallel computers allows large-scale inverse problems, comprising several hundred thousand nodes with hundreds of sources and receivers, to be solved. The new method is tested using computer-generated data from two- and three-dimensional synthetic models. For each inversion a choice of penalty parameters, gauging the level of model covariance information imposed, has to be made and the level of regularization required is hard to estimate. We find that running a suite of inversions with varying penalty parameters and subsequent examination of the results (including inspection of residual maps) offers a viable method for choosing appropriate numerical values for the penalty levels. In the applications we found the inversion process to be highly non-linear. Inversion models from intermediate steps of the iterative inversion show structure in places that do not exhibit structure in the true model and only at later iterations do anomalies move to the correct location in the modelling domain. This result indicates that linearized inversions that fail to re-linearize during the inversion process will fail to find meaningful inversion images. The inversion images achieved using the new method recover the important features of the true models, including the approximate magnitudes of the conductivity anomalies and the magnitudes and directions of anisotropy anomalies. The inversion images are generally 'blurred', that is sharp edges are smoothed, and the recovered magnitudes of conductivity, anisotropy and anisotropy direction are generally under-estimated.