We have developed a new uncertainty estimation method that accounts for nonlinearity inherent in most geophysical problems, allows for the explicit search of model posterior space, is scalable, and maintains computational efficiencies on the order of deterministic inverse solutions. We accomplish this by combining an efficient parameter reduction technique, a parameter constraint mapping routine, a sparse geometric sampling scheme, and an efficient forward solver. In order to reduce our model domain and determine an independent basis, we implement both a typical principal component analysis, which factorizes the model covariance matrix, and an alternative compression method, based on singular-value decomposition, which acts on training models, directly, and is storage efficient. Once we have a reduced base, we map parameter constraints, from our original model domain, to this reduced domain to define a bounded geometric region of feasible model space. We utilize an optimal scheme to sample this reduced model space that uses Smolyak sparse grids and minimizes the number of forward solves by finding the sparsest sampling required to produce convergent uncertainty measures. The result is an ensemble of equivalent models, consistent with our inverse solution structure, which is used to infer inverse uncertainty. We tested our method with a 1D synthetic example, a comparison with a published Metropolis-Hastings sampling example, and an extension to 2D problems using a field data inversion.