BUG: MLLE implementation does not yield expected results

### Describe the bug

The implemented manifold learning method Modified Locally Linear Embedding (MLLE) in /sklearn/manifold/_locally_linear.py line 403 ('method==modified') does not produce the expected results, see attached img.

From the reference [3] it is clear that MLLE should reconstruct a lower dimesional manifold which preserves local geometric characteristics of the higher-dimensional input data.
This means data points of the same color which are close in the 'Original S-curve' 3D data should also be close in reconstructed 2D manifold. 
Also points which are far away from each other along the local coordinate along the S-shape  should be far away in the reconstruction. Meaning the yellow points and the blue points should be the furthest apart in the reconstructed data.

As can be seen from the [example](https://scikit-learn.org/stable/auto_examples/manifold/plot_compare_methods.html) the current implementation of MLLE does not statisfy these properties.

[3] Zhang, Z. & Wang, J. MLLE: Modified Locally Linear Embedding Using Multiple Weights. (http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.70.382)
![MLLE_Bugreport](https://user-images.githubusercontent.com/61986870/188318658-081ce71d-6e11-47a1-a158-c67f07adc632.jpg)


### Steps/Code to Reproduce

To reproduce, run the corresponding part of the example [Comparison of Manifold Learning methods](https://scikit-learn.org/stable/auto_examples/manifold/plot_compare_methods.html):

```
from numpy.random import RandomState
import matplotlib.pyplot as plt
from matplotlib import ticker

# unused but required import for doing 3d projections with matplotlib < 3.2
import mpl_toolkits.mplot3d

from sklearn import manifold, datasets

rng = RandomState(0)

n_samples = 1500
S_points, S_color = datasets.make_s_curve(n_samples, random_state=rng)

def plot_3d(points, points_color, title):
    x, y, z = points.T

    fig, ax = plt.subplots(
        figsize=(6, 6),
        facecolor="white",
        tight_layout=True,
        subplot_kw={"projection": "3d"},
    )
    fig.suptitle(title, size=16)
    col = ax.scatter(x, y, z, c=points_color, s=50, alpha=0.8)
    ax.view_init(azim=-60, elev=9)
    ax.xaxis.set_major_locator(ticker.MultipleLocator(1))
    ax.yaxis.set_major_locator(ticker.MultipleLocator(1))
    ax.zaxis.set_major_locator(ticker.MultipleLocator(1))

    fig.colorbar(col, ax=ax, orientation="horizontal", shrink=0.6, aspect=60, pad=0.01)
    plt.show()


def plot_2d(points, points_color, title):
    fig, ax = plt.subplots(figsize=(3, 3), facecolor="white", constrained_layout=True)
    fig.suptitle(title, size=16)
    add_2d_scatter(ax, points, points_color)
    plt.show()


def add_2d_scatter(ax, points, points_color, title=None):
    x, y = points.T
    ax.scatter(x, y, c=points_color, s=50, alpha=0.8)
    ax.set_title(title)
    ax.xaxis.set_major_formatter(ticker.NullFormatter())
    ax.yaxis.set_major_formatter(ticker.NullFormatter())


plot_3d(S_points, S_color, "Original S-curve samples")

n_neighbors = 12  # neighborhood which is used to recover the locally linear structure
n_components = 2  # number of coordinates for the manifold

params = {
    "n_neighbors": n_neighbors,
    "n_components": n_components,
    "eigen_solver": "auto",
    "random_state": rng,
}

lle_standard = manifold.LocallyLinearEmbedding(method="standard", **params)
S_standard = lle_standard.fit_transform(S_points)

lle_ltsa = manifold.LocallyLinearEmbedding(method="ltsa", **params)
S_ltsa = lle_ltsa.fit_transform(S_points)

lle_hessian = manifold.LocallyLinearEmbedding(method="hessian", **params)
S_hessian = lle_hessian.fit_transform(S_points)

lle_mod = manifold.LocallyLinearEmbedding(method="modified", modified_tol=0.8, **params)
S_mod = lle_mod.fit_transform(S_points)

fig, axs = plt.subplots(
    nrows=2, ncols=2, figsize=(7, 7), facecolor="white", constrained_layout=True
)
fig.suptitle("Locally Linear Embeddings", size=16)

lle_methods = [
    ("Standard locally linear embedding", S_standard),
    ("Local tangent space alignment", S_ltsa),
    ("Hessian eigenmap", S_hessian),
    ("Modified locally linear embedding", S_mod),
]
for ax, method in zip(axs.flat, lle_methods):
    name, points = method
    add_2d_scatter(ax, points, S_color, name)

plt.show()
```

### Expected Results

The expected result is a manifold which preserves local geometric features, similar to other methods such as  'Isomap', 'LTSA' or 'Hessian Eigenmap'

The results shown here are generated with my own implementation of the paper by Zhang, Z. & Wang, J. 'MLLE: Modified Locally Linear Embedding Using Multiple Weights'. However, the notation is from their later publication [5], where they call the method 'NEML', but it is the same algorithm.

To implement replace the code of the current MLLE implementation with the following snippet:
```
from scipy.linalg import schur

elif method == "modified":
        reg = 1e-3
        roh = np.zeros((N))
        eigenValues = []
        eigenVectors = []
        WOpti = []

        if n_neighbors < n_components:
            raise ValueError("method neml requires n_neighbors >= n_components")

            # step 1
        J = nbrs.kneighbors(
            X, n_neighbors=n_neighbors + 1, return_distance=False
        )

        for i in range(N):
            # step 1.1
            xi = X[i, :]
            Ji = J[i, 1:]  # neighbrhood of xi without xi itself.

            # step 1.2
            Gi = X[Ji, :] - xi
            C = Gi @ Gi.transpose()
            C_tilde = C + np.eye(n_neighbors, n_neighbors) * reg * np.trace(C)
            wiOpt = np.linalg.solve(C_tilde, np.ones((n_neighbors, 1)))
            wiOpt = wiOpt / np.sum(wiOpt)  # regularization as suggested at p.3

            # step 1.3
            [S, V] = schur(C, output="real")
            ei = np.sort(np.diag(S))
            ei = ei[::-1]
            ei[ei <= np.finfo(float).eps] = 0
            JIi = np.argsort(np.diag(S))
            JIi = JIi[::-1]
            roh[i] = np.sum(ei[n_components : n_neighbors]) / np.sum(ei[0 : n_components])

            # save eigenvalues and vectors for next steps
            eigenValues.append(ei)
            eigenVectors.append(V[:, JIi])
            WOpti.append(wiOpt)

        # step 2
        rohSorted = np.sort(roh)
        eta = rohSorted[int(np.ceil(N / 2))]
        s = np.zeros(N)
        for i in range(N):
            l = n_neighbors - n_components
            lambdas = eigenValues[i]
            while (
                np.sum(lambdas[n_neighbors - l :]) / np.sum(lambdas[: n_neighbors - l]) > eta
                and l > 1
            ):
                l = l - 1

            s[i] = l

        # step 3
        M = np.zeros((N, N))
        for i in range(N):
            Ji = J[i, 1:]
            Vi = eigenVectors[i]
            Vhat = Vi[:, int(n_neighbors - s[i]) :]
            vi = np.sum(Vhat, 0)  # equivalent to Vhat.transpose() @ np.ones(k, 1)
            alphai = np.linalg.norm(vi) / np.sqrt(s[i])
            ui = np.ones((int(s[i]))) * alphai - vi
            uiNorm = np.linalg.norm(ui)
            if uiNorm > 1e-5:
                ui = ui / uiNorm
            else:
                ui = np.zeros((int(s[i]), 1))
            Hi = np.eye(int(s[i])) - 2 * np.outer(ui, ui)
            Wi = (1 - alphai) ** 2 * WOpti[i] @ np.ones((1, int(s[i]))) + (
                2 - alphai
            ) * Vhat @ Hi

            # build Alignment Matrix
            M[i, i] = M[i, i] + s[i]
            M[np.ix_(Ji, Ji)] = M[np.ix_(Ji, Ji)] + Wi @ Wi.transpose()
            M[Ji, i] = M[Ji, i] - np.sum(Wi, 1)
            M[i, Ji] = M[Ji, i]

        if M_sparse:
            M = csr_matrix(M)
```
[5] Wang J., Z. & Zhang, Z. Nonlinear embedding preserving muliple local-linearities. Pattern Recognition.  43 (2010)

![MLLE_ExpectedResult](https://user-images.githubusercontent.com/61986870/188319603-b65f6826-c789-4659-91ac-71a96fff3234.PNG)


### Actual Results

The acutal result is a construction which does not seem to preserve the desired local gemeometric characteristics.

![MLLE_Result](https://user-images.githubusercontent.com/61986870/188319560-a3a6426a-8200-4bff-af19-88574fe7619c.PNG)


### Versions

```shell
System:
    python: 3.8.10 (default, Jun 22 2022, 20:18:18)  [GCC 9.4.0]
executable: /home/markuswnuk/projects/scikit-learn/sklearn-env/bin/python3
   machine: Linux-5.10.102.1-microsoft-standard-WSL2-x86_64-with-glibc2.29

Python dependencies:
      sklearn: 1.2.dev0
          pip: 20.0.2
   setuptools: 65.3.0
        numpy: 1.23.2
        scipy: 1.9.1
       Cython: 0.29.32
       pandas: 1.4.4
   matplotlib: 3.5.3
       joblib: 1.1.0
threadpoolctl: 3.1.0

Built with OpenMP: True

threadpoolctl info:
       user_api: blas
   internal_api: openblas
         prefix: libopenblas
       filepath: /home/markuswnuk/projects/scikit-learn/sklearn-env/lib/python3.8/site-packages/numpy.libs/libopenblas64_p-r0-742d56dc.3.20.so
        version: 0.3.20
threading_layer: pthreads
   architecture: Haswell
    num_threads: 8

       user_api: blas
   internal_api: openblas
         prefix: libopenblas
       filepath: /home/markuswnuk/projects/scikit-learn/sklearn-env/lib/python3.8/site-packages/scipy.libs/libopenblasp-r0-9f9f5dbc.3.18.so
        version: 0.3.18
threading_layer: pthreads
   architecture: Haswell
    num_threads: 8

       user_api: openmp
   internal_api: openmp
         prefix: libgomp
       filepath: /usr/lib/x86_64-linux-gnu/libgomp.so.1.0.0
        version: None
    num_threads: 8
```

COMMENT:
Thank you for pointing this out @markuswnuk91. I think this is just due to the use of `modified_tol=0.8` in the plotting demo. Allowing it to use the default value of `1e-12` instead of manually setting it to ` 0.8` seems to fix this:

![download](https://user-images.githubusercontent.com/34613774/188617474-33cc96fd-1e83-4a50-a475-fd5893d90715.png)


Would you be interested in opening a PR to correct the plotting demo?

Labels: Bug, module:manifold