Indoor positioning systems cannot rely on conventional localization methods, such as GPS, to locate devices because of interference with the structure of buildings. One solution is to use magnetic positioning, which is based on spatial variations in the patterns of the ambient magnetic field. To model magnetic fields, Gaussian process regression is used, providing predictions of the magnetic field at unvisited locations along with uncertainty quantification. These predictions and their uncertainties are valuable information for probabilistic localization algorithms used for magnetic positioning. Full Gaussian process regression has poor scalability, becoming computationally intractable from roughly 10,000 one-dimensional measurements due to its associated cubic computational complexity. In the existing literature, approximations for Gaussian process regression have been extensively studied to reduce this computational complexity. Of these approximations, only approximations involving basis functions and local experts have been used in the context of scalable magnetic field modeling. A favorable approximation framework from existing literature uses structured kernel interpolation (SKI), allowing for fast regression through efficient Krylov subspace methods. The SKI framework is favorable as it allows for fast regression in low dimensions without introducing boundary effects. In this thesis, the SKI framework is used to approximate two distinct magnetic field models: the shared model, which considers independence between the magnetic field components with shared hyperparameters, and the scalar potential model, which includes physical properties of the magnetic field (Maxwell’s equations) in the model. The scalability of the approach is shown using simulations and experiments with magnetic field measurements. Through the simulations, it is shown the SKI framework accurately and efficiently approximates the models. The applicability of the SKI framework for scalable magnetic field modeling is investigated using data collected using a motion capture suit, the Xsens MVN Link Suit. In the final experiment, a magnetic field map is constructed based on more than 40,000 three-dimensional measurements without splitting the data set, which took less than one minute on a standard laptop.