Large numbers of high-dimensional patterns are collected in a variety of disciplines, from astronomy to bioinformatics. In this article, we present an approach to non-linear dimensionality reduction based on fitting nearest neighbor regression to the unsupervised regression framework for learning of low-dimensional manifolds. For each high-dimensional pattern, a low-dimensional latent point is generated. The dimensionality of the induced optimization problem grows with the number of patterns. To cope with the large solution space, an iterative solution construction scheme is proposed. In this paper, we introduce two strategies to embed high-dimensional data. First, the latent sorting approach allows embeddings in a one-dimensional latent space corresponding to a sorting of the high-dimensional patterns. Second, Gaussian embeddings randomly generate candidate positions based on sampling from the Gaussian distribution employing distances on data space as variances. Kernel functions increase the flexibility of the approach by mapping the patterns to feature spaces. We analyze and compare the algorithms experimentally on a set of test functions.