We present a method for inverting surface magnetic data to recover 3-D susceptibility models. To allow the maximum flexibility for the model to represent geologically realistic structures, we discretize the 3-D model region into a set of rectangular cells, each having a constant susceptibility. The number of cells is generally far greater than the number of the data available, and thus we solve an underdetermined problem. Solutions are obtained by minimizing a global objective function composed of the model objective function and data misfit. The algorithm can incorporate a priori information into the model objective function by using one or more appropriate weighting functions. The model for inversion can be either susceptibility or its logarithm. If susceptibility is chosen, a positivity constraint is imposed to reduce the nonuniqueness and to maintain physical realizability. Our algorithm assumes that there is no remanent magnetization and that the magnetic data are produced by induced magnetization only. All minimizations are carried out with a subspace approach where only a small number of search vectors is used at each iteration. This obviates the need to solve a large system of equations directly, and hence earth models with many cells can be solved on a deskside workstation. The algorithm is tested on synthetic examples and on a field data set.