Our study presents a probabilistic non-rigid point set registration method to deal with large and uneven deformations. Our method treats the registration as a density estimation problem. In our method, we add two key constraints to enforce landmark correspondences and preserve local neighborhood structure. We assume that the landmarks, which represent the salient points in the point sets, are given or can be detected using keypoint detectors such as scale-invariant feature transform or MeshDOG. By enforcing landmark correspondences, we preserve the overall global shape of the point set with significant deformations. Furthermore, by leveraging stochastic neighbor embedding, we incorporate constraints to preserve local neighborhood structure, which penalizes incoherent transformation within a neighborhood. Our experimental results in both 2D and 3D datasets show that our method outperforms state-of-the-art methods in a large degree of deformations. In particular, quantitative results show that the error is 29% better than the second-best result (from the state-of-the-art methods). Our analysis shows that a relatively small number of landmarks is sufficient to deal with large deformations. Finally, our study shows that our method is computationally comparable to state-of-the-art methods.