Manifold Learning with Isomap¶
Important: Please read the installation page for details about how to install the toolboxes. $\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{\umin}[1]{\underset{#1}{\min}\;}$ $\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{\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.
from __future__ import division
import numpy as np
import scipy as scp
import pylab as pyl
import matplotlib.pyplot as plt
from nt_toolbox.general import *
from nt_toolbox.signal import *
import warnings
warnings.filterwarnings('ignore')
%matplotlib inline
%load_ext autoreload
%autoreload 2
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*(.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.
D1 = np.repeat(np.sum(X**2, 0)[:,np.newaxis], n, 1)
D1 = np.sqrt(D1 + np.transpose(D1) - 2*np.dot(np.transpose(X), X))
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()