```python
import numpy as np
from sklearn.manifold import SpectralEmbedding
from sklearn.datasets import make_swiss_roll, make_s_curve

import matplotlib.pyplot as plt
import mpl_toolkits.mplot3d.axes3d as p3
plt.style.use('ggplot')

%matplotlib inline
```

This notebook will walk through the steps of Laplacian Eigenmaps (LE) algorithm. It will be a step-by-step walkthrough of the algorithm and towards the notebook, I will work on some common known speed-up attempts.


```python
seed = 123
rng = np.random.seed(123)
n_samples = 1500
noise = 0.1
random_state = seed

data, color = make_swiss_roll(n_samples=n_samples, noise=noise, random_state=random_state)
data, color = make_s_curve(n_samples=n_samples, noise=noise, random_state=random_state)
```

\begin{align}
    g &= \int_a^b f(x)dx \label{eq1} \\
    a &= b + c \label{eq2}
\end{align}

See (\ref{eq1})


```python
fig = plt.figure()

ax = p3.Axes3D(fig)

ax.scatter(data[:,0], data[:, 1], data[:,2], c=color, cmap=plt.cm.Spectral)
ax.set_title("Original Data")
plt.show()
```

### Laplacian Eigenmaps (Sklearn)


```python
%%time

n_components = 2
affinity = 'nearest_neighbors'
n_neighbors = 10
n_jobs = -1
random_state = 123
eigen_solver = 'arpack'

# initialize le model
le_model = SpectralEmbedding(
    n_components=n_components,
    affinity=affinity,
    n_neighbors=n_neighbors,
    n_jobs=n_jobs,
    random_state=random_state,
    eigen_solver=eigen_solver
)

# fit and transform data
embedding = le_model.fit_transform(data);
```

    CPU times: user 175 ms, sys: 15.9 ms, total: 191 ms
    Wall time: 214 ms



```python
fig, ax = plt.subplots(figsize=(10, 6))

ax.spy(le_model.affinity_matrix_, markersize=1.0)
plt.show()
```


```python
print(embedding.shape)
```

    (1500, 2)



```python
fig, ax = plt.subplots()

ax.scatter(embedding[:, 0], embedding[:, 1], c=color)
ax.set_title('Projected Data')
plt.show()
```

## Adjacency Matrix Construction

### Nearest Neighbours Search


```python
# some baseline parameters
n_neighbors = 10
algorithm = 'brute'
metric = 'euclidean'
p=2
n_jobs = -1

# initialize nn model
nn_model = NearestNeighbors(
    n_neighbors=n_neighbors,
    metric=metric,
    algorithm=algorithm,
    p=p,
    n_jobs=n_jobs
)

# fit nn model to data
nn_model.fit(data);

# grab distances and indices
dists, indices = nn_model.kneighbors(
    data, 
    n_neighbors=n_neighbors, 
    return_distance=True
)
```

### Weighted Distances


```python
# Heat kernel
def heat_kernel(distances, length_scale=None):
    if length_scale is None:
        length_scale = 1.0
#         length_scale = np.sqrt(distances.shape[1] / 2.0)
#         length_scale = 1.0 / distances.shape[1]
#     return np.exp(- length_scale * distances**2)
    return np.exp(- distances**2 / length_scale)
```


```python
# transform distances with heat kernel
dists = heat_kernel(dists)
```

### Construct Graph


```python
# Construct sparse KNN Graph
n_samples = data.shape[0]

indptr = np.arange(0, n_samples * n_neighbors + 1, n_neighbors)

adjacency_matrix = csr_matrix((dists.ravel(), indices.ravel(), indptr), shape=(n_samples, n_samples))

# ensure that its symmetrix
adjacency_matrix = 0.5 * (adjacency_matrix + adjacency_matrix.T)
```

#### Peak at Adjacency Matrix


```python

fig, ax = plt.subplots()

ax.spy(adjacency_matrix, markersize=1.0)
ax.set_title('Adjacency Matrix', pad=15.0)

plt.show()
```

## Laplacian Matrix

Some notes about some different Laplacian matrices:

**Unnormalized Graph Laplacian**

$$L=D-W$$

**Symmetric Normalized Graph Laplacian**

$$L_{Sym}=D^{-\frac{1}{2}}LD^{-\frac{1}{2}}$$

**Random Walk Normalized Laplacian**

$$L_{rw}=D^{-1}L$$

**Random Walk Transition Matrix**

$$L_{rwt}=D^{-1}W$$

### Laplacian and Degree Matrix


```python
def graph_laplacian(adjacency_matrix, graph_type='normalized', return_diag=True):
    
    n_samples = adjacency_matrix.shape[0]
    
    # Get degree vector
    degree = np.array(adjacency_matrix.sum(axis=1)).squeeze()
    
    # Create Sparse matrix for degree
    degree_mat = spdiags(
        degree, diags=0, m=n_samples, n=n_samples
    )
    
    if graph_type in ['unnormalized']:
        laplacian = degree_mat - adjacency_matrix
    elif graph_type in ['normalized']:
        laplacian = degree_mat - adjacency_matrix
        norm = 1 / np.sqrt(degree_mat)
        laplacian = norm @ adjacency_matrix @ norm
    

    
    # Set all values 
    
    # set diagonal elements of array to negative
    
    
    return laplacian, degree_mat
```


```python
%%time
laplacian, degree = graph_laplacian(adjacency_matrix)
```


```python
# print(degree)
```


```python

```


```python
print(degree.shape)
print(laplacian.diagonal().shape)
```

    (1500,)
    (1500,)



```python
laplacian, degree = graph_laplacian(adjacency_matrix, normed=True, return_diag=True)
```


```python
print(laplacian.shape, degree.shape)
```

    (1500, 1500) (1500,)



```python
fig, ax = plt.subplots(nrows=1, ncols=2)

ax[0].spy(laplacian, markersize=1.0)
ax[0].set_title('Laplacian Matrix', pad=15.0)

ax[1].spy(degree, markersize=1.0)
ax[1].set_title('Degree Matrix', pad=15.0)

plt.show()
```

## EigenValue Decomposition

|        Algorithm       | Laplacian Equation      | Generalized Eigenvalue | Standard Eigenvalue |
|:----------------------:|-------------------------|:----------------------:|---------------------|
| Unnormalized Laplacian |        $$L=D-W$$        |          $$$$          |                     |
| Normalized Laplacian   | $$L=D^{-1/2}LD^{-1/2}$$ |          2.75          |                     |
| Random Walk            |                         |          2.03          |                     |
| ReNormalized           |                         |          0.64          |                     |
| Geometric              |                         |                        |                     |

### Generalized Eigenvalue Solver


```python
%%time

# Transform equation
laplacian_sol = -1 * laplacian

n_components = 2
solver = 'LM'      # Smallest to Largest
sigma = 1.0
eigen_tol = 0.0
v0 = np.random.uniform(-1, 1, laplacian.shape[0])

eigenvalues, eigenvectors = eigsh(
    laplacian_sol,
    k=n_components,
    which=solver,
    sigma=sigma,
    tol=eigen_tol,
    v0=v0
)

# Transform eigenvectors
embedding = eigenvectors.T[n_components::-1]

```

    CPU times: user 141 ms, sys: 3.27 ms, total: 144 ms
    Wall time: 83.6 ms



```python
print(eigenvalues.shape, eigenvectors.shape)
```

    (2,) (1500, 2)



```python
fig, ax = plt.subplots()

ax.scatter(embedding[0, :], embedding[1,:], c=color)
ax.set_title('Projected Data')
plt.show()
```


```python

```

1cfa3cf2b8ba00333fa6560d4ed46ce04dbe0e3e
ipynb
Jupyter Notebook
notebooks/uncategorized/eigenmap/le_sklearn.ipynb
J-Garcke-SCAI/jaxkern
9de7ebf52fe2d186d316350a6692b2ecc0885adc
MIT
2020-09-28T07:39:16.000Z
2022-03-11T14:09:41.000Z
notebooks/uncategorized/eigenmap/le_sklearn.ipynb
J-Garcke-SCAI/jaxkern
9de7ebf52fe2d186d316350a6692b2ecc0885adc
MIT
2020-09-25T01:25:57.000Z
2020-10-09T16:15:49.000Z
notebooks/uncategorized/eigenmap/le_sklearn.ipynb
J-Garcke-SCAI/jaxkern
9de7ebf52fe2d186d316350a6692b2ecc0885adc
MIT
2021-05-25T21:59:58.000Z
2022-01-11T07:23:32.000Z
Qwen/Qwen-72B
1. YES
2. YES