3. Differential Forms

Differential forms provide a coordinate-free framework for integration on manifolds. A differential form has a degree $k$, corresponding to the dimension of the objects it can be integrated over. In particular, a 1-form integrates along curves, assigning a real number to each oriented path, while a 2-form integrates over surfaces, assigning a real number to each oriented surface. More generally, a $k$-form can be integrated over $k$-dimensional submanifolds.

Given a manifold $M$, a $k$-form $\omega$ is an alternating $k$-linear map $$\omega_p:(T_p M)^k\to \mathbb{R},$$ for each point $p$ on $M$. This is an element of $\Lambda^kT_p^*M$.

At each point $p\in M$, the metric $g_p:T_pM \times T_pM\to \mathbb{R}$ induces an identification $$\flat :T_pM\to T^*_pM, \ \ \ v\mapsto g_p(v,\cdot). $$ The inverse is denoted by $\#$ and the pair of isomorphisms $(\flat, \#)$ are termed to be the musical isomorphisms. The Riemannian metric extends to an inner product on exterior powers: $$g^{(k)}:\Lambda^kT_pM\times \Lambda^kT_pM\to\mathbb{R} $$ defined on decomposable elements by $$g^{(k)}\!\left(v_1 \wedge \cdots \wedge v_k,\;w_1 \wedge \cdots \wedge w_k\right)=\det\!\bigl( g(v_i, w_j) \bigr).$$

This induces a musical isomorphism $$\flat : \Lambda^k T_p M \longrightarrow \Lambda^k T_p^* M,\qquad\alpha \longmapsto g^{(k)}(\alpha, \cdot).$$ Hence, a $k$-form may be equivalently regarded as a metric-dual of a $k$-vector field.

We will exactly use this duality to define differential forms. On the generators, we have that $\#: f \nabla h \mapsto f dh$, with the defining relation that on the ambient coordinate function we have $dx_j(x_i) = \delta_{ij}$. As such, we use exactly the same representaion for 1-forms as for vector fields: $$\alpha = \sum_{i=1}^{n_1}\sum^d_{j=1}\alpha_{ij} \phi_i dx_j.$$

We write $k$-forms as sums $$\beta = \sum_i f_i dh^1_i\wedge \dots \wedge dh_i^k.$$ Using the ambient coordinate functions, we get the expression $$\beta = \sum_{1\leq j_1<\dots <j_j\leq d} f_{j_1,\dots,j_k} dx_{j_1} \wedge \dots\wedge dx_{j_k}=\sum_{J=1}^{\binom{d}{k}} f_J dx_J, $$ where $J=(j_1,\dots, j_k)$ is a multi-index in lexicographic (increasing) order.

Riemannian Metric

The metric on $k$-forms uses the determinant to compute a signed volume $$g(fdh_1\wedge\dots\wedge dh_k, f' dh_1'\wedge\dots\wedge dh'_k)=ff'\mathrm{det}(\Gamma(h_i,h_j'))$$

Visualising Differential Forms

Since 1-forms are dual to vector fields, we visualise the corresponding vector field.

To compute these objects, we generalise the 1-form approach by evaluating the $k$-form $\alpha$ against coordinate vector fields. This produces a $d \times \dots \times d$ $k$-tensor that is skew-symmetric in its indices:

$$g(\alpha, dx^{j_1} \wedge \dots \wedge dx^{j_k})_{j_1,\dots,j_k=1,\dots,d} \in \mathbb{R}^{d \times \dots \times d}$$

2-forms in 2D

A 2-form in 2d gives a skew symmetric $2\times 2$ matrix $$\begin{pmatrix} 0 & a \\ -a & 0 \end{pmatrix},$$ where $a\in\mathbb{R}$ is the signed area magnitude, which we can visualise as a coloured disk of size $|a|$ and signed colour. The following code block constructs the first few generating 2-forms on a point clouds, and visualises them with the above prescription.

The DiffusionGeometry backend class handles the creation of $k$-forms. The data can be either provided in the point wise ambient basis, or in the spectral coefficients basis. When it comes to the visualisation, the `visualisation' submodule contains all the required functions. The following is a minimal example.

import os, sys
sys.path.insert(0, os.path.abspath('..'))  # Add parent directory to path

from diffusion_geometry import DiffusionGeometry
from diffusion_geometry.visualisation import clean_fig, plot_2form_2d

import numpy as np

from plotly.subplots import make_subplots



# Visualise the first 9 basis 2-forms on a grid

n = 200
# Create a meshgrid in 2D
x = np.linspace(-1, 1, int(np.sqrt(n)))
y = np.linspace(-1, 1, int(np.sqrt(n)))
X, Y = np.meshgrid(x, y) 
data = np.vstack([X.ravel(), Y.ravel()]).T

dg = DiffusionGeometry.from_point_cloud(data)

specs = [[{"type": "xy"} for _ in range(3)] for _ in range(3)]
fig = make_subplots(
    rows=3,
    cols=3,
    specs=specs,
    horizontal_spacing=0.05,
    vertical_spacing=0.05,
)

for i in range(3):
    for j in range(3):
        # Creates a 2-form with all coefficients zero
        two_form = dg.form_space(2).zeros()

        # Set the coefficient of the (i,j)-th basis 2-form to 1
        basis_index = i * 3 + j  # Assuming a 3x3 grid of basis forms
        two_form.coeffs[basis_index] = 1

        # Plot the 2-form on the disc
        fig_2d = plot_2form_2d(data, two_form.to_ambient())

        for trace in fig_2d.data:
            fig.add_trace(trace, row=i+1, col=j+1)
        

fig.update_layout(
    width=600,
    height=600,
    margin=dict(l=0, r=0, t=0, b=0),
)
clean_fig(fig)
fig.show()

2-forms in 3D

These yield a $3 \times 3$ skew-symmetric matrix. By orthogonally diagonalising the matrix, we find the specific 2D plane where the form acts: $$\begin{pmatrix} 0 & a & 0 \\ -a & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}$$ Here, $a$ is the signed area magnitude within the plane spanned by the first two eigenvectors. This is all handled by the function 'plot_2form_3d' from the visualisation module.

3-forms in 3D

These yield a $3 \times 3 \times 3$ skew-symmetric tensor. It can be represented as a collection of matrices where $a$ is the signed volume magnitude: $$\left[ \begin{pmatrix} 0 & 0 & 0 \\ 0 & 0 & a \\ 0 & -a & 0 \end{pmatrix}, \begin{pmatrix} 0 & 0 & -a \\ 0 & 0 & 0 \\ a & 0 & 0 \end{pmatrix}, \begin{pmatrix} 0 & a & 0 \\ -a & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix} \right]$$ These computations are handled by the function 'plot_3form_3d' from the visualisation module.

from plotly.subplots import make_subplots

from diffusion_geometry.visualisation import plot_3form_3d, plot_scatter_3d
from figures.generate_data import gen_3d_data


n = 1000
data = gen_3d_data(kind = 'ball', minor_radius = 1.0, major_radius=2.0, n = n, noise = 0.0, seed = 0)[0]
data = data[:, [1, 2, 0]]
dg = DiffusionGeometry.from_point_cloud(data)

num_rows = 2
num_cols = 4

specs = [[{"type": "scene"} for _ in range(num_cols)] for _ in range(num_rows)]
fig = make_subplots(
    rows=num_rows,
    cols=num_cols,
    specs=specs,
    horizontal_spacing=0.0,
    vertical_spacing=0.0,
)


n_sphere = 12
x_vals_all, y_vals_all, z_vals_all = [], [], []

camera = dict(eye=dict(x=-1., y=-1.2, z=0.3),
                center=dict(x=0, y=0, z=0),
                up=dict(x=0, y=0, z=1))


# ---------- build subplots ----------
for col_idx in range(1, num_cols + 1):
    # function example
    f = dg.function_space.zeros()
    f.coeffs[col_idx - 1] = 1
    fig1 = plot_scatter_3d(data, color=f.to_ambient(), size=2)

    # 3-form example
    a = dg.form_space(3).zeros()
    a.coeffs[col_idx - 1] = 1
    fig2 = plot_3form_3d(data, a.to_ambient(), base_size=0.0, size_scale=30.0, camera=camera)

    # add to subplot grid
    fig.add_traces(list(fig1.data), rows=[1]*len(fig1.data), cols=[col_idx]*len(fig1.data))
    fig.add_traces(list(fig2.data), rows=[2]*len(fig2.data), cols=[col_idx]*len(fig2.data))

    # collect coordinates from both figures
    for subfig in (fig1, fig2):
        for t in subfig.data:
            if hasattr(t, "x") and t.x is not None:
                x_vals_all.extend(t.x)
            if hasattr(t, "y") and t.y is not None:
                y_vals_all.extend(t.y)
            if hasattr(t, "z") and t.z is not None:
                z_vals_all.extend(t.z)

# ---------- compute shared ranges ----------
xrange_3d = [min(x_vals_all), max(x_vals_all)]
yrange_3d = [min(y_vals_all), max(y_vals_all)]
zrange_3d = [min(z_vals_all), max(z_vals_all)]

# ---------- synchronize all 3D axes ----------
camera = dict(
    center=dict(x=0, y=0, z=0),
    eye=dict(x=-1., y=-1.1, z=0.3),
    up=dict(x=0, y=0, z=1),
)

clean_fig(fig)

# apply same ranges and camera to every scene in layout
for scene_key in [k for k in fig.layout if k.startswith("scene")]:
    fig.layout[scene_key].update(
        aspectmode="data",
        camera=camera,
        xaxis=dict(range=xrange_3d),
        yaxis=dict(range=yrange_3d),
        zaxis=dict(range=zrange_3d),
    )

# ---------- layout and export ----------
fig.update_layout(
    width=1000,
    height=500,
    margin=dict(l=0, r=0, t=0, b=0),
)
fig.show()

Wedge Product

Differential forms have a notion of multiplication that extends the usual multiplication of functions (i.e. 0-forms).

If $\alpha \in \Omega^k(M)$ and $\beta \in \Omega^l(M)$, we can form their wedge product $\alpha \wedge \beta \in \Omega^{k+l}(M)$ by

$$(f\, dh_1 \wedge \dots \wedge dh_k) \wedge (f' \, dh_1' \wedge \dots \wedge dh_l') = ff' \, dh_1 \wedge \dots \wedge dh_k \wedge dh_1' \wedge \dots \wedge dh_l'.$$

In our discretisation, we can compute the wedge product of

$$\boldsymbol{\alpha} = \sum_{i=1}^{n_1} \sum_{J = 1}^{\binom{d}{k}} \boldsymbol{\alpha}_{iJ} \, p_i \, dx_{(J)}, \qquad\boldsymbol{\beta} = \sum_{i'=1}^{n_1} \sum_{J' = 1}^{\binom{d}{l}} \boldsymbol{\beta}_{i'J'} \, p_{i'} \, dx_{(J')}$$

as

$$\boldsymbol{\alpha} \wedge \boldsymbol{\beta}= \sum_{i,i'=1}^{n_1} \sum_{J = 1}^{\binom{d}{k}}\sum_{J' = 1}^{\binom{d}{l}} \boldsymbol{\alpha}_{iJ} \, \boldsymbol{\beta}_{i'J'} \, p_i \, p_{i'} \, dx_{(J)} \wedge dx_{(J')}.$$

Note that the multi-indices $J = (j_1,\dots,j_k)$ and $J' = (j'_1,\dots,j'_l)$ are in lexicographic (increasing) order. However, while it is true that $dx_{(J)} \wedge dx_{(J')} = dx_{(j_1,...,j_k,j'_1,...,j'_l)}$, the concatenated multi-index $(j_1,...,j_k,j'_1,...,j'_l)$ is not necessarily in lexicographic order. We therefore need to sort this list and multiply the result by the sign of the permutation.

The class 'Form' handles this, and we can wedge two forms by using the ^ operator.

# Generate figure visualising the wedge product of two 1-forms in 2D


from diffusion_geometry.tensors.functions.function import Function
from diffusion_geometry.visualisation import plot_quiver_2d


def genereate_1form_1form_figure(fig, current_row):
    r=1
    n_side = 28
    x = np.linspace(-r, r, n_side)
    y = np.linspace(-r, r, n_side)
    X, Y = np.meshgrid(x, y)
    data = np.column_stack((X.ravel(), Y.ravel())) * [1.3, 1]
    # data += np.random.randn(data.shape[0], 2) * 0.01
   

    dg = DiffusionGeometry.from_point_cloud(data)

    # Define two scalar functions that are cosines and sines
    def f(x,y):
        return np.sin(5*x)*np.cos(4*x)
    def g(x,y):
        return np.cos(3*y)*np.sin(2*y)


    f_data = f(data[:,0], data[:,1])
    g_data = g(data[:,0], data[:,1])

    f = Function.from_pointwise_basis(f_data, dg)
    g = Function.from_pointwise_basis(g_data, dg)

    # Construct 1-forms by taking the exterior derivative
    alpha = f.d()
    beta = g.d()

    # Compute the wedge product. Note that this can simply be written as alpha ^ beta since we have overloaded the ^ operator
    wedged_2_form = alpha ^ beta

    # Plot settings
    quiver_scale = 0.05
    arrow_scale = 0.5
    line_width = 1

    # Plot the three components: alpha, beta, and alpha ^ beta
    fig1 = plot_quiver_2d(data, alpha.to_ambient(), scale=quiver_scale, arrow_scale=arrow_scale, line_width=line_width)
    fig2 = plot_quiver_2d(data, beta.to_ambient(), scale=quiver_scale, arrow_scale=arrow_scale, line_width=line_width)
    fig3 = plot_2form_2d(data, wedged_2_form.to_ambient(), n_circle=dg.n)

    fig.add_traces(list(fig1.data), rows=current_row, cols=1)
    fig.add_traces(list(fig2.data), rows=current_row, cols=2)
    fig.add_traces(list(fig3.data), rows=current_row, cols=3)

fig = make_subplots(
    rows=1,
    cols=3,
    specs=[[{"type": "xy"} for _ in range(3)] for _ in range(1)],
    horizontal_spacing=0.05,
    vertical_spacing=0.1,
)

genereate_1form_1form_figure(fig, current_row=1)

clean_fig(fig)
fig.update_layout(
    width=800,
    height=400, 
)




fig.update_yaxes(scaleanchor="x", scaleratio=1)
fig.update_xaxes(scaleanchor="y", scaleratio=1)
fig.show()