The construction of smooth potential field models has many geophysical applications. The recently-developed method of harmonic splines produces magnetic field models at the core surface which are maximally smooth in the sense of minimization of certain special norms. They do not exhibit the highly oscillatory fields produced by models derived from a least squares analysis with a truncated spherical harmonic series. Modeling the data by harmonic splines requires solving a square system of equations with dimension equal to the number of data. Too many data have been collected since the 1960s for this method to be practical. We produce almost optimally smooth models by the following method. Since each spline function for the optimal model corresponds to an observation location (called a knot), we select a subset of these splines with knots well-distributed around the earth's surface. In this depleted basis we then find the smoothest model subject to an appropriate fit to all of the data. This reduces the computational problem to one comparable to least-squares analysis while nearly preserving the optimality inherent in the original harmonic spline models.