It might be possible to solve your example problem using a procedure similar to nonclassical metric MDS (using the stress criterion). Initialize the 'projected points' to lie on a sphere (more on this later). Then, use an optimization solver to find the projected points that minimize the objective function. There are a few differences compared to ordinary MDS. 1) Geodesic distances must be computed between projected points, rather than Euclidean distances. This is easy when the manifold is a sphere. 2) The optimization must obey the constraint that the projected points lie on a sphere.
Fortunately, there are existing, open source solvers for performing optimization where the parameters are constrained to lie on a particular manifold. Since the parameters here are the projected points themselves, this is exactly what's needed. Manopt is a package for Matlab, and Pymanopt is for Python. They include spheres as one of the supported manifolds, and others are available.
The quality of the final result will depend on the initialization. This is also the case for ordinary, nonclassical MDS, where a good initial configuration is often obtained using classical MDS (which can be solved efficiently as an eigenvalue problem). For 'spherical MDS', you could take the following approach for initialization. Perform ordinary MDS, isomap, or some other nonlinear dimensionality reduction technique to obtain coordinates in a Euclidean space. Then, map the resulting points onto the surface of a sphere using a suitable projection. For example, to project onto a 3-sphere, first perform ordinary dimensionality reduction to 2d. Map the resulting points onto a 3-sphere using something like a stereographic projection. If the original data lies on some manifold that's topologically equivalent to a sphere, then it might be more appropriate to perform initial dimensionality reduction to 3d (or do nothing if they're already in 3d), then normalize the vectors to pull them onto a sphere. Finally, run the optimization. As with ordinary, nonclassical MDS, multiple runs can be performed using different initial conditions, then the best result selected.
It should be possible to generalize to other manifolds, and to other objective functions. For example, we could imagine converting the objective functions of other nonlinear dimensionality reduction algorithms to work on spheres or other manifolds.