Laplacians of simplicial complexes

In [1]:
%matplotlib inline
In [2]:
import numpy as np
from scipy import sparse
import scipy.sparse.linalg
import matplotlib as mpl
from matplotlib import pyplot as plt
import gudhi
import pygsp as pg

import sys
sys.path.append('..')
from data.s2_6_complex_to_laplacians import build_boundaries, build_laplacians

1 Build an alpha complex from a point cloud

The simplicial complex is built from a Delaunay triangulation. See the CGAL documentation.

Alternatives:

In [3]:
n_points = 100  # Approximate.
dim = 2

rs = np.random.RandomState(None)
points = np.concatenate([
    [0.2, 1] * rs.uniform(size=(n_points//4, dim)),
    [0.2, 0.8] + [0.6, 0.2] * rs.uniform(size=(n_points//6, dim)),
    [0.2, 0] + [0.6, 0.2] * rs.uniform(size=(n_points//6, dim)),
    [0.8, 0] + [0.2, 1] * rs.uniform(size=(n_points//4, dim)),
    [0.4, 0.4] + [0.2, 0.2] * rs.uniform(size=(n_points//6, dim))
])

points = np.random.uniform(size=(n_points, dim))
#points = pg.graphs.Grid2d(10).coords

n_points = points.shape[0]
print(f'{n_points} points')

plt.scatter(*points.T);
100 points

from scipy import spatial

tri = spatial.Delaunay(points) print(f'{tri.simplices.shape[0]} triangles')

plt.triplot(points[:,0], points[:,1], tri.simplices);

In [4]:
ac = gudhi.AlphaComplex(points)
st = ac.create_simplex_tree()

before = st.num_simplices()
_ = st.prune_above_filtration(1e-2)
print(f'filtration: {before} => {st.num_simplices()} simplices')

assert st.num_vertices() == n_points
assert st.dimension() == dim
filtration: 573 => 498 simplices

2 Extract simplices

  • $n_0$ is the number of 0-simplices (nodes)
  • $n_1$ is the number of 1-simplices (edges)
  • $n_2$ is the number of 2-simplices (triangles)
  • $n_k$ is the number of $k$-simplices

simplices is a list of dictionaries, with one dictionary per degree $k$. Each dictionary maps a simplex (represented as a set of vertices) to an integer that will be its index in the boundary and Laplacian operators.

In [5]:
def extract_simplices(simplex_tree):
    simplices = [dict() for _ in range(simplex_tree.dimension()+1)]
    for simplex, _ in simplex_tree.get_skeleton(simplex_tree.dimension()):
        k = len(simplex)
        simplices[k-1][frozenset(simplex)] = len(simplices[k-1])
    return simplices

simplices = extract_simplices(st)
In [6]:
for k, s in enumerate(simplices):
    print(f'n_{k} = {len(s):,} {k}-simplices')
n_0 = 100 0-simplices
n_1 = 251 1-simplices
n_2 = 147 2-simplices

2.1 Plotting

In [7]:
def get_positions(simplices, dim):
    polygons = list()
    for i, simplex in enumerate(simplices[dim].keys()):
        assert simplices[dim][simplex] == i  # Dictionary is ordered.
        polygon = list()
        for vertex in simplex:
            polygon.append(points[vertex])
        polygons.append(polygon)
    return polygons

lines = get_positions(simplices, 1)
triangles = get_positions(simplices, 2)
In [8]:
def value2color(values):
    values -= values.min()
    values /= values.max()
    return mpl.cm.viridis(values)

def plot_nodes(colors, ax=None, **kwargs):
    if ax is None:
        fig, ax = plt.subplots()
    ax.scatter(points[:, 0], points[:, 1], c=colors, **kwargs)
    return ax.figure, ax

def plot_edges(colors, ax=None, **kwargs):
    if ax is None:
        fig, ax = plt.subfigs()
    colors = value2color(colors)
    collection = mpl.collections.LineCollection(lines, colors=colors, **kwargs)
    ax.add_collection(collection)
    ax.autoscale()

def plot_triangles(colors, ax=None, **kwargs):
    if ax is None:
        fig, ax = plt.subfigs()
    colors = value2color(colors)
    for triangle, color in zip(triangles, colors):
        triangle = plt.Polygon(triangle, color=color, **kwargs)
        ax.add_patch(triangle)
    ax.autoscale()
    
def plot_triangles_plain( ax=None, **kwargs):
    if ax is None:
        fig, ax = plt.subfigs()
    colors = len(simplices[1])*['linen']
    for triangle, color in zip(triangles, colors):
        triangle = plt.Polygon(triangle, color=color, **kwargs)
        ax.add_patch(triangle)
    ax.autoscale()

3 Build boundary operators

$$ \nabla_k \in \mathbb{R}^{n_{k-1} \times n_k}, \ k \in [1, K] $$
  • Also known as incidence matrices.
  • Each $k$-simplex is connected to its $(k-1)$-faces.
  • They are the gradient operators (and their transpose are the divergence operators).
In [9]:
boundaries = build_boundaries(simplices)

3.1 Discrete derivatives and integrals

  • $\nabla_k^\top \in \mathbb{R}^{n_k \times n_{k-1}}$ is the difference operator (discrete exterior derivative)
  • $\langle \nabla_k^\top f, g \rangle = (\nabla_k^\top f)^\top g = f^\top \nabla_k g = \langle f, \nabla_k g \rangle$, hence $\nabla_k$ is the codifferential of $\nabla_k^\top$. The codifferential is the adjoint of the difference.
  • The codifferential represents sums ("finite integrals") over $k$-surfaces.
In [10]:
s1 = np.zeros(len(simplices[1]))
s1[100:150] = 1

#s1 = np.random.uniform(size=len(simplices[1]))

s0 = boundaries[0] @ s1
s2 = boundaries[1].T @ s1

fig, ax = plt.subplots()
plot_nodes(s0, ax, zorder=3)
plot_edges(s1, ax, zorder=2)
plot_triangles(s2, ax, zorder=1)

#plt.colorbar(ax.collections[0], ax=ax)

4 Build Laplacians

$$ L_k = L_k^\text{down} + L_k^\text{up} = \nabla_{k}^\top \nabla_{k} + \nabla_{k+1} \nabla_{k+1}^\top \in \mathbb{R}^{n_k \times n_k}, \ k \in [0, K] $$
  • Graph Laplacian $L_0 = \nabla_{1} \nabla_{1}^\top$
  • Last Laplacian $L_K = \nabla_{K}^\top \nabla_{K}$
In [11]:
laplacians = build_laplacians(boundaries)

for k, laplacian in enumerate(laplacians):
    print('{}-simplices: {:,} simplices, {:.2%} sparse'.format(k, laplacian.shape[0], laplacian.nnz/np.prod(laplacian.shape)))
    assert laplacian.shape == (len(simplices[k]), len(simplices[k]))
0-simplices: 100 simplices, 6.02% sparse
1-simplices: 251 simplices, 2.59% sparse
2-simplices: 147 simplices, 2.49% sparse

5 Eigenvectors

$$ L_k = V_k \Lambda_k V_k^\top, \ V_k = [v_1, \dots, v_{n_k}] \in \mathbb{R}^{n_k \times n_k} $$

5.1 Vertices (0-simplices)

In [12]:
def plot_eigen(laplacian, plot, n_eigenvectors=4):
    values, vectors = sparse.linalg.eigsh(laplacian, n_eigenvectors, which='SM')
    fix, axes = plt.subplots(1, n_eigenvectors+3, figsize=(20, 3))
    axes[0].plot(values, '.')
    axes[1].scatter(vectors[:, 0], vectors[:, 1])
    axes[2].scatter(vectors[:, 2], vectors[:, 3])
    for i in range(n_eigenvectors):
        plot(vectors[:, i], axes[i+3])
        axes[i+3].axis('off')

plot_eigen(laplacians[0], plot_nodes)

5.2 Edges (1-simplices)

The edges indeed seem to be embeded in linear subspaces.

TODO:

In [13]:
plot_eigen(laplacians[1], plot_edges)

5.3 Triangles (2-simplices)

In [14]:
plot_eigen(laplacians[2], plot_triangles)

5.4 Heat diffusion

We study the following dynamical system.

$$\frac{\partial s(t)}{\partial t}=-L_k s(t)$$
In [15]:
values, vectors = [], []
for j , laplacian in enumerate(laplacians):
    w, v = np.linalg.eigh(laplacians[j].toarray())
    values.append(w)
    vectors.append(v)
In [ ]:
 
In [16]:
class Heat:

    def __init__(self, t, lmax):
        self.t = t
        self.lmax = lmax

    def __call__(self, x):
        return np.exp(-self.t * x / self.lmax)

    def plot(self):
        x = np.linspace(0, self.lmax)
        plt.plot(x, self(x))
dim=0
heat = Heat(t=10, lmax=values[dim][-1])
heat.plot()
In [17]:
s = np.zeros(len(simplices[0]))
s[50:70] = 1

#s = np.random.uniform(size=len(simplices[0]))

def filter(d, signal):
    return vectors[d] @ (heat(values[d]) * (vectors[d].T @ signal))

s = filter(0, s)
print(s.shape)

fig, ax = plot_nodes(s)
plt.colorbar(ax.collections[0], ax=ax)
#plot_edges(np.zeros(len(simplices[1])), ax=ax, zorder=0)

print(np.linalg.norm(s, 1))
(100,)
20.000000620898522

Heat diffusion on 1-simplices using the 1-Laplacian

In [18]:
dim=1
iterations=1000
heat = Heat(t=iterations, lmax=values[dim][-1])
heat.plot()
In [19]:
s = np.zeros(len(simplices[dim]))
#s[100:150] =1
s = np.random.uniform(size=len(simplices[dim]))
init_cond=np.copy(s)
def filter(d, signal):
    return vectors[d] @ (heat(values[d]) * (vectors[d].T @ signal))

s = filter(dim, s)
filtred_s=np.copy(s)
print(s.shape)
fig, ax = plt.subplots()
plot_edges(s, ax,norm=False, zorder=2)
plot_triangles_plain(ax, zorder=1)
#plt.colorbar(ax.collections[0], ax=ax)
print(np.linalg.norm(s, 1))
(251,)
139.46565826729983

Convergence of the dynamical system.

As studied in Control Using Higher Order Laplacians in Network Topologies, Chapter IV.

  • In case $H_1=0$ the dynamical system converges to 0.
  • In case $H_1\neq 0$ the dynamical system converges to some point in $ker(L_1)$
In [20]:
if np.isclose(values[1][0],10**-10):
    epsilon=0.0001
    harm_dist=np.linalg.norm(vectors[1][0]-filtred_s, 2)
    c=~np.isclose(values[1],10**-8)
    t_epsilon=-1/values[1][c][0]*np.log(epsilon/(np.linalg.norm(init_cond, 2)*len(init_cond)))
    print('The solution is epsilon-harmonic for all t > {}'.format(t_epsilon))
    print('After {} iterations the distance from one harmonic representative is {}'.format(iterations,harm_dist))
else:
    print('The solution converges to zero. After {} iterations the norm is {}'.format(iterations,np.linalg.norm(filtred_s, 2)))
The solution is epsilon-harmonic for all t > 183.92344235057902
After 1000 iterations the distance from one harmonic representative is 1.368369778518699
In [21]:
graph = pg.graphs.Grid2d(10)
g = pg.filters.Heat(graph, 100)
s = np.zeros(graph.N)
s[50:70] = 1
s = g.filter(s, method='exact')
graph.plot_signal(s)
print(np.linalg.norm(s, 1))
2020-11-04 19:48:30,224:[WARNING](pygsp.graphs.graph._check_fourier_properties): The eigenvalues vector G.e is not available, we need to compute the Fourier basis. Explicitly call G.compute_fourier_basis() once beforehand to suppress the warning.
19.999999999999993

6 Test

Simplicial complex and ground truths from Control Using Higher Order Laplacians in Network Topologies, figure 3.

In [22]:
simplices = [
    {
        frozenset([0]): 0,
        frozenset([1]): 1,
        frozenset([2]): 2,
        frozenset([3]): 3,
        frozenset([4]): 4,
    },
    {
        frozenset([0, 1]): 0,
        frozenset([0, 2]): 1,
        frozenset([1, 2]): 2,
        frozenset([1, 3]): 3,
        frozenset([2, 3]): 4,
        frozenset([2, 4]): 5,
        frozenset([3, 4]): 6,
    },
    {
        frozenset([1, 2, 3]): 0,
    },
]

boundaries = build_boundaries(simplices)
laplacians = build_laplacians(boundaries)

_, vectors = sparse.linalg.eigsh(laplacians[1], 2, which='SM')
print(vectors)

values, vectors = np.linalg.eigh(laplacians[1].toarray())
plt.plot(values, '.');
[[ 0.4999086   0.36199972]
 [-0.49990845 -0.36199963]
 [ 0.432163    0.06092054]
 [ 0.06774557  0.30107912]
 [-0.36441737  0.24015856]
 [ 0.2966719  -0.54123765]
 [-0.29667178  0.5412376 ]]