Manifold Learning with Isomap¶
Recover a low-dimensional representation of a curved point cloud using neighborhood-graph distances. Compare geodesic and Euclidean distances, then use a spectral embedding to reveal the manifold’s intrinsic geometry.
Run this tour¶
Run the cells in order with a Python 3 kernel. The first cell locates the companion data and toolbox and installs missing dependencies when needed. All worked examples include their implementation directly in this notebook. Random seeds make comparisons reproducible; you can change them to explore other samples.
# Locate the companion toolbox locally, or fetch it for a standalone/Colab copy.
from pathlib import Path
import importlib.util
import os
import subprocess
import sys
working = Path.cwd()
candidates = [working, working / "python", working.parent / "python"]
python_dir = next((p for p in candidates if (p / "nt_toolbox").is_dir()), None)
if python_dir is None:
checkout = working / "numerical-tours-support"
if not checkout.exists():
subprocess.run(
[
"git",
"clone",
"--depth",
"1",
"--branch",
"master",
"https://github.com/gpeyre/numerical-tours.git",
str(checkout),
],
check=True,
)
python_dir = checkout / "python"
os.chdir(python_dir)
if str(python_dir) not in sys.path:
sys.path.insert(0, str(python_dir))
requirements = python_dir / "requirements.txt"
if any(
importlib.util.find_spec(name) is None
for name in [
"numpy",
"scipy",
"matplotlib",
"skimage",
"sklearn",
"pywt",
"ipywidgets",
"cvxpy",
"skfmm",
"autograd",
"progressbar",
"celer",
]
):
subprocess.run(
[sys.executable, "-m", "pip", "install", "-r", str(requirements)], check=True
)
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(0)
plt.rcParams.update(
{
"figure.figsize": (8, 4),
"figure.dpi": 100,
"axes.spines.top": False,
"axes.spines.right": False,
"font.size": 11,
"image.cmap": "gray",
}
)
%matplotlib inline
$\newcommand{\dotp}[2]{\langle #1, #2 \rangle}$ $\newcommand{\enscond}[2]{\lbrace #1, #2 \rbrace}$ $\newcommand{\pd}[2]{ \frac{ \partial #1}{\partial #2} }$ $\newcommand{\umin}[1]{\underset{#1}{\min}\;}$ $\newcommand{\umax}[1]{\underset{#1}{\max}\;}$ $\newcommand{\uargmin}[1]{\underset{#1}{argmin}\;}$ $\newcommand{\norm}[1]{\|#1\|}$ $\newcommand{\abs}[1]{\left|#1\right|}$ $\newcommand{\choice}[1]{ \left\{ \begin{array}{l} #1 \end{array} \right. }$ $\newcommand{\pa}[1]{\left(#1\right)}$ $\newcommand{\diag}[1]{{diag}\left( #1 \right)}$ $\newcommand{\qandq}{\quad\text{and}\quad}$ $\newcommand{\qwhereq}{\quad\text{where}\quad}$ $\newcommand{\qifq}{ \quad \text{if} \quad }$ $\newcommand{\qarrq}{ \quad \Longrightarrow \quad }$ $\newcommand{\ZZ}{\mathbb{Z}}$ $\newcommand{\CC}{\mathbb{C}}$ $\newcommand{\RR}{\mathbb{R}}$ $\newcommand{\EE}{\mathbb{E}}$ $\newcommand{\Zz}{\mathcal{Z}}$ $\newcommand{\Ww}{\mathcal{W}}$ $\newcommand{\Vv}{\mathcal{V}}$ $\newcommand{\Nn}{\mathcal{N}}$ $\newcommand{\NN}{\mathcal{N}}$ $\newcommand{\Hh}{\mathcal{H}}$ $\newcommand{\Bb}{\mathcal{B}}$ $\newcommand{\Ee}{\mathcal{E}}$ $\newcommand{\Cc}{\mathcal{C}}$ $\newcommand{\Gg}{\mathcal{G}}$ $\newcommand{\Ss}{\mathcal{S}}$ $\newcommand{\Pp}{\mathcal{P}}$ $\newcommand{\Ff}{\mathcal{F}}$ $\newcommand{\Xx}{\mathcal{X}}$ $\newcommand{\Mm}{\mathcal{M}}$ $\newcommand{\Ii}{\mathcal{I}}$ $\newcommand{\Dd}{\mathcal{D}}$ $\newcommand{\Ll}{\mathcal{L}}$ $\newcommand{\Tt}{\mathcal{T}}$ $\newcommand{\si}{\sigma}$ $\newcommand{\al}{\alpha}$ $\newcommand{\la}{\lambda}$ $\newcommand{\ga}{\gamma}$ $\newcommand{\Ga}{\Gamma}$ $\newcommand{\La}{\Lambda}$ $\newcommand{\Si}{\Sigma}$ $\newcommand{\be}{\beta}$ $\newcommand{\de}{\delta}$ $\newcommand{\De}{\Delta}$ $\newcommand{\phi}{\varphi}$ $\newcommand{\th}{\theta}$ $\newcommand{\om}{\omega}$ $\newcommand{\Om}{\Omega}$
This tour explores the Isomap algorithm for manifold learning.
The <http://waldron.stanford.edu/~isomap/ Isomap> algorithm is introduced in
A Global Geometric Framework for Nonlinear Dimensionality Reduction, J. B. Tenenbaum, V. de Silva and J. C. Langford, Science 290 (5500): 2319-2323, 22 December 2000.
import numpy as np
import scipy as scp
import pylab as pyl
import matplotlib.pyplot as plt
from nt_toolbox.general import np, plt, pylab, rescale
from nt_toolbox.signal import np, plt, pylab
import warnings
%matplotlib inline
Graph Approximation of Manifolds¶
Manifold learning consist in approximating the parameterization of a manifold represented as a point cloud.
First we load a simple 3D point cloud, the famous Swiss Roll.
Number of points.
n = 1000
Random position on the parameteric domain.
from numpy import random
x = random.rand(2, n)
Mapping on the manifold.
v = 3 * np.pi / 2 * (0.1 + 2 * x[0, :])
X = np.zeros([3, n])
X[1, :] = 20 * x[1, :]
X[0, :] = -np.cos(v) * v
X[2, :] = np.sin(v) * v
Parameter for display.
ms = 200
el = 20
az = -110
Display the point cloud.
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure(figsize=(15, 11))
ax = fig.add_subplot(111, projection="3d")
# swiss roll
ax.scatter(
X[0, :],
X[1, :],
X[2, :],
c=plt.cm.jet((X[0, :] ** 2 + X[2, :] ** 2) / 100),
s=ms,
lw=0,
alpha=1,
)
# params
ax.set_xlim(np.min(X[0, :]), np.max(X[0, :]))
ax.set_ylim(np.min(X[1, :]), np.max(X[1, :]))
ax.set_zlim(np.min(X[2, :]), np.max(X[2, :]))
ax.axis("off")
ax.view_init(elev=el, azim=az)
Compute the pairwise Euclidean distance matrix.
from scipy.spatial.distance import cdist
D1 = cdist(X.T, X.T)
Number of NN for the graph.
k = 6
Compute the k-NN connectivity.
DNN, NN = np.sort(D1), np.argsort(D1)
NN = NN[:, 1 : k + 1]
DNN = DNN[:, 1 : k + 1]
Adjacency matrix, and weighted adjacency.
from scipy import sparse
B = np.tile(np.arange(0, n), (k, 1))
A = sparse.coo_matrix((np.ones(k * n), (np.ravel(B, order="F"), np.ravel(NN))))
Weighted adjacency (the metric on the graph).
W = sparse.coo_matrix((np.ravel(DNN), (np.ravel(B, order="F"), np.ravel(NN))))
Display the graph.
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure(figsize=(15, 11))
ax = fig.add_subplot(111, projection="3d")
# swiss roll
ax.scatter(
X[0, :],
X[1, :],
X[2, :],
c=plt.cm.jet((X[0, :] ** 2 + X[2, :] ** 2) / 100),
s=ms,
lw=0,
alpha=1,
)
# graph
I, J, V = sparse.find(A)
xx = np.vstack((X[0, I], X[0, J]))
yy = np.vstack((X[1, I], X[1, J]))
zz = np.vstack((X[2, I], X[2, J]))
for i in range(len(I)):
ax.plot(xx[:, i], yy[:, i], zz[:, i], color="black")
# params
ax.axis("off")
ax.set_xlim(np.min(X[0, :]), np.max(X[0, :]))
ax.set_ylim(np.min(X[1, :]), np.max(X[1, :]))
ax.set_zlim(np.min(X[2, :]), np.max(X[2, :]))
ax.view_init(elev=el, azim=az)
plt.show()
Floyd Algorithm to Compute Pairwise Geodesic Distances¶
A simple algorithm to compute the geodesic distances between all pairs of points on a graph is Floyd iterative algorithm. Its complexity is $\mathcal O(n^3)$ where $n$ is the number of points. It is thus quite slow for sparse graph, where Dijkstra runs in $\mathcal O(n^2\log(n))$.
Floyd algorithm iterates the following update rule, for $k=1,\dots,n$
$D(i,j) \leftarrow \min(D(i,j), D(i,k)+D(k,j))$,
with the initialization $D(i,j)=W(i,j)$ if $W(i,j)>0$, and $D(i,j)=Inf$ if $W(i,j)=0$.
Make the graph symmetric.
# An edge is present when either point includes the other among its neighbors.
W = W.tocsr().maximum(W.T.tocsr())
D = W.toarray()
Initialize the matrix.
D[D == 0] = float("inf")
Add connexion between a point and itself.
np.fill_diagonal(D, 0)
Worked example 1
Implement the Floyd algorithm to compute the full distance matrix $D$, where $D(i,j)$ is the geodesic distance between
from scipy.sparse.csgraph import shortest_path
D = shortest_path(W, directed=False)
if not np.isfinite(D).all():
raise ValueError("The neighborhood graph is disconnected; increase k.")
assert np.allclose(D, D.T)
Find index of vertices that are not connected to the main manifold.
Iremove = np.flatnonzero(~np.isfinite(D[:, 0]))
Remove Inf remaining values (disconnected components).
assert np.isfinite(D).all(), "Every pair needs a finite geodesic distance."
Isomap with Classical Multidimensional Scaling¶
Isomap perform the dimensionality reduction by applying multidimensional scaling.
Please refers to the tours on Bending Invariant for detail on Classical MDS (strain minimization).
Worked example 2
Perform classical MDS to compute the 2D flattening.
from numpy import linalg
J = np.identity(n) - np.ones([n, n]) / n
K = -1 / 2.0 * np.dot(J, np.dot(D**2, J))
[val, Xstrain] = sparse.linalg.eigsh(K, k=2, which="LA")
Xstrain = Xstrain * np.tile(np.sqrt(np.maximum(val, 0)), (n, 1))
Xstrain = np.real(np.transpose(Xstrain))
# plot size
plt.figure(figsize=(15, 6))
# plot points
plt.scatter(
Xstrain[0, :],
Xstrain[1, :],
ms,
c=plt.cm.jet((X[0, :] ** 2 + X[2, :] ** 2) / 100),
lw=0,
alpha=1,
)
# plot vertices
I, J, V = sparse.find(A)
xx = np.vstack((Xstrain[0, I], Xstrain[0, J]))
yy = np.vstack((Xstrain[1, I], Xstrain[1, J]))
for i in range(len(I)):
plt.plot(xx[:, i], yy[:, i], color="black")
# params
plt.axis("off")
plt.xlim(np.min(Xstrain[0, :] - 1), np.max(Xstrain[0, :]) + 1)
plt.ylim(np.min(Xstrain[1, :] - 1), np.max(Xstrain[1, :]) + 1)
plt.show()
Redess the points using the two leading eigenvectors of the covariance matrix (PCA correction).
[L, U] = linalg.eigh(np.dot(Xstrain, np.transpose(Xstrain)) / n)
Xstrain1 = np.dot(np.transpose(U), Xstrain)
Remove problematic points.
Xstrain1[:, Iremove] = float("inf")
Display the final result of the dimensionality reduction.
# plot size
plt.figure(figsize=(15, 6))
# plot points
plt.scatter(
Xstrain1[0, :],
Xstrain1[1, :],
ms,
c=plt.cm.jet((X[0, :] ** 2 + X[2, :] ** 2) / 100),
lw=0,
alpha=1,
)
# plot vertices
I, J, V = sparse.find(A)
xx = np.vstack((Xstrain1[0, I], Xstrain1[0, J]))
yy = np.vstack((Xstrain1[1, I], Xstrain1[1, J]))
for i in range(len(I)):
plt.plot(xx[:, i], yy[:, i], color="black")
# params
plt.axis("off")
plt.xlim(np.min(Xstrain1[0, :] - 1), np.max(Xstrain1[0, :]) + 1)
plt.ylim(np.min(Xstrain1[1, :] - 1), np.max(Xstrain1[1, :]) + 1)
plt.show()
For comparison, the ideal locations on the parameter domain.
Y = np.vstack((v, X[1, :]))
Y[0, :] = rescale(Y[0, :], min(Xstrain[0, :]), max(Xstrain[0, :]))
Y[1, :] = rescale(Y[1, :], min(Xstrain[1, :]), max(Xstrain[1, :]))
Display the ideal graph on the reduced parameter domain.
# plot size
plt.figure(figsize=(15, 6))
# plot points
plt.scatter(
Y[0, :],
Y[1, :],
ms,
c=plt.cm.jet((X[0, :] ** 2 + X[2, :] ** 2) / 100),
lw=0,
alpha=1,
)
# plot vertices
I, J, V = sparse.find(A)
xx = np.vstack((Y[0, I], Y[0, J]))
yy = np.vstack((Y[1, I], Y[1, J]))
for i in range(len(I)):
plt.plot(xx[:, i], yy[:, i], color="black")
# params
plt.axis("off")
plt.xlim(np.min(Y[0, :] - 1), np.max(Y[0, :]) + 1)
plt.ylim(np.min(Y[1, :] - 1), np.max(Y[1, :]) + 1)
plt.show()
Isomap with SMACOF Multidimensional Scaling¶
It is possible to use SMACOF instead of classical scaling.
Please refers to the tours on Bending Invariant for detail on both Classical MDS (strain minimization) and SMACOF MDS (stress minimization).
Worked example 3
Perform stress minimization MDS using SMACOF to compute the 2D flattening.
fig = plt.figure(figsize=(15, 15))
niter = 150
stress = []
Xstress = np.copy(X)
ndisp = [1, 5, 10, min(niter, 100), np.inf]
k = 0
for i in range(niter):
if ndisp[k] == i:
ax = fig.add_subplot(2, 2, k + 1, projection="3d")
# swiss roll
ax.scatter(
Xstress[0, :],
Xstress[1, :],
Xstress[2, :],
c=plt.cm.jet((X[0, :] ** 2 + X[2, :] ** 2) / 100),
s=ms,
lw=0,
alpha=1,
)
# graph
I, J, V = sparse.find(A)
xx = np.vstack((Xstress[0, I], Xstress[0, J]))
yy = np.vstack((Xstress[1, I], Xstress[1, J]))
zz = np.vstack((Xstress[2, I], Xstress[2, J]))
for i in range(len(I)):
ax.plot(xx[:, i], yy[:, i], zz[:, i], color="black")
# params
ax.axis("off")
ax.set_xlim(np.min(Xstress[0, :]), np.max(Xstress[0, :]))
ax.set_ylim(np.min(Xstress[1, :]), np.max(Xstress[1, :]))
ax.set_zlim(np.min(Xstress[2, :]), np.max(Xstress[2, :]))
ax.view_init(elev=el, azim=az)
k += 1
# Compute the distance matrix.
D1 = np.repeat(np.sum(Xstress**2, 0)[:, np.newaxis], n, 1)
D1 = np.sqrt(
np.maximum(
D1 + np.transpose(D1) - 2 * np.dot(np.transpose(Xstress), Xstress), 0
)
)
# Compute the scaling matrix.
B = -D / np.maximum(D1, 1e-10 * np.ones(np.shape(D1)))
B = B - np.diag(np.sum(B, 0))
# update
Xstress = np.transpose(np.dot(B, np.transpose(Xstress))) / n
# Xstress = Xstress-repmat(mean(Xstress,2), [1 n]);
# record stress
stress = stress + [np.sqrt(np.sum(abs(D - D1) ** 2) / n**2)]
plt.show()
Plot stress evolution during minimization.
plt.figure(figsize=(10, 7))
plt.plot(stress, ".-")
plt.show()
Compute the main direction of the point clouds.
[L, U] = linalg.eigh(np.dot(Xstress, np.transpose(Xstress)) / n)
[L, I] = np.sort(L), np.argsort(L)
U = U[:, I[1:3]]
Project the points on the two leading eigenvectors of the covariance matrix (PCA final projection).
Xstress1 = np.dot(np.transpose(U), Xstress)
Remove problematic points.
Xstress1[:, Iremove] = np.inf
Display the final result of the dimensionality reduction.
# plot size
plt.figure(figsize=(15, 6))
# plot points
plt.scatter(
Xstress1[1, :],
Xstress1[0, :],
ms,
c=plt.cm.jet((X[0, :] ** 2 + X[2, :] ** 2) / 100),
lw=0,
alpha=1,
)
# plot vertices
I, J, V = sparse.find(A)
xx = np.vstack((Xstress1[1, I], Xstress1[1, J]))
yy = np.vstack((Xstress1[0, I], Xstress1[0, J]))
for i in range(len(I)):
plt.plot(xx[:, i], yy[:, i], color="black")
# params
plt.axis("off")
plt.xlim(np.min(Xstress1[1, :] - 1), np.max(Xstress1[1, :]) + 1)
plt.ylim(np.min(Xstress1[0, :] - 1), np.max(Xstress1[0, :]) + 1)
plt.show()
References and further reading¶
James A. Sethian. A Fast Marching Level Set Method for Monotonically Advancing Fronts. 1996, PNAS 93(4), 1591–1595. A causal numerical solver for the eikonal equation.
Edsger W. Dijkstra. A Note on Two Problems in Connexion with Graphs. 1959, Numerische Mathematik 1, 269–271. Shortest paths with nonnegative edge weights.
Stanley Osher and James A. Sethian. Fronts Propagating with Curvature-Dependent Speed: Algorithms Based on Hamilton–Jacobi Formulations. 1988, Journal of Computational Physics 79(1), 12–49. Implicit front propagation using level-set functions.
Joshua B. Tenenbaum, Vin de Silva, and John C. Langford. A Global Geometric Framework for Nonlinear Dimensionality Reduction. 2000, Science 290, 2319–2323. Geodesic distances and multidimensional scaling for manifold learning.
Gabriel Peyré. Advanced Signal, Image and Surface Processing. 2010, course notes. A mathematical companion to the Numerical Tours.