We consider the problem of nonlinear dimensionality reduction: given a training set of high-dimensional data whose "intrinsic" low dimension is assumed known, find a feature extraction map to low-dimensional space, a reconstruction map back to high-dimensional space, and a geometric description of the dimension-reduced data as a smooth manifold. We introduce a complexity-regularized quantization approach for fitting a Gaussian mixture model to the training set via a Lloyd algorithm. Complexity regularization controls the trade-off between adaptation to the local shape of the underlying manifold and global geometric consistency. The resulting mixture model is used to design the feature extraction and reconstruction maps and to define a Riemannian metric on the low-dimensional data. We also sketch a proof of consistency of our scheme for the purposes of estimating the unknown underlying pdf of high-dimensional data.