2. Vector Fields

The Carré du Champ Identity

The fundamental object of Riemannian geometry is the Riemannian metric. Given two vector fields $X=\sum_{i=1}^d f_i \nabla x_i$ and $Y=\sum_{j=1}^d h_j \nabla x_j$ in $\mathbb{R}^d$, the Riemannian metric of $X$ and $Y$ decomposes as $$g(X,Y) = \sum_{i,j=1}^df_i h_j g(\nabla x_i, \nabla x_j).$$ The terms $g(\nabla x_i, \nabla x_j)$ are called the carré du champ of $x_i$ and $x_j$, and denoted $\Gamma(x_i, x_j)$. Crucially, these terms can be expressed in terms of the diffusion process. As diffusion generalises to point clouds, we are able to extend the computation of the carré du champ from manifolds to point clouds. It can be computed in a data-driven way via the Markov chain $P_{ij}$ defined earlier: $$\Gamma_i(f,h) = \frac{1}{2\rho(x_i)^2} \sum_{j=1}^n P_{ij} \left(f_j - \sum^n_{k=1}P_{jk}f_k\right) \left(h_j - \sum^n_{k=1}P_{jk}h_k\right).$$

Inspired by this, we define vector fields to be sums $$X= \sum_{j=1}^d f_j \nabla x_j.$$ We now decompose the function $f_j$, however, for computational reason we introduce another truncation parameter $n_1\leq n_0$ that is able to compress the function space for vector fields even more. We have: $$X = \sum_{i=1}^{n_1}\sum_{j=1}^d X_{ij} \phi_i \nabla x_j.$$ In this way, vector fields can be represented by the $n_1\times d$ coefficient matrices $(X_{ij})_{i\leq n_1, j\leq d}$.

Riemannian Metric

We can define a pointwise inner product between two vector fields $X$ and $Y$ by setting $$g(f\nabla h, f' \nabla h'),$$ and extending by linearity to arbitrary vector fields.

To visualise a vector field $X$, we compute its projections onto the ambient spanning set of vector fields in $\mathbb{R}^d$. This means computing $$(g(X, \nabla x_1), \dots, g(X, \nabla x_d))\in \mathbb{R}^d. $$

Completely analagous to function, we can define vector fields pointwise as $n\times d$ matrices. The framework then handles the conversion into the spectral basis behind the scenes. In the code below, we define a vector field in the ambient space that points along the $x$-axis.

import sys
import os
sys.path.append(os.path.abspath('..'))  # Add parent directory to path

from diffusion_geometry import DiffusionGeometry
from diffusion_geometry import VectorField

import numpy as np
from diffusion_geometry.visualisation import *
from plotly.subplots import make_subplots


np.random.seed(0)

# Generate our running example: synthetic 2D point cloud M
n = 150
xx = np.linspace(0, 2*np.pi, n)
M = np.array([[np.cos(x), np.sin(x)] for x in xx]) + 0.2*np.random.randn(xx.shape[0], 2)
M = 1*np.array([[np.cos(x)/(1 + np.sin(x)**2), np.sin(x)*np.cos(x)/(1 + np.sin(x)**2)] for x in xx]) + 0.1*M + 0.0*np.random.randn(xx.shape[0], 2) 
M = np.concatenate((M, 0.18*np.random.randn(200,2) - [0.2,0]))

# Truncate function space to first 50 eigenfunction 
# and truncate vector field coefficients to first 40 eigenfunction
dg = DiffusionGeometry.from_point_cloud(M, n_function_basis=50, n_coefficients=40)

# Define a vector field in pointwise basis: constant field in x-direction
X_pointwise = np.zeros_like(M)
X_pointwise[:, 0] = 1

# Construct VectorField object from pointwise basis representation
X = VectorField.from_pointwise_basis(X_pointwise, dg)

# Visualise the vector field, careful to convert to ambient coordinates
scale = 0.1
line_width = 1
plot_quiver_2d(dg.immersion_coords, X.to_ambient(), scale=scale, line_width=line_width).show()

Both the classes 'Function' and 'VectorField' inherit from a parent class 'Tensor'. The following example shows how to take the Riemannian metric between two vector fields. Inside each tangent space, this gives us the projection of one vector field onto the other. On a disc, we now define a rotational and an outward pointing radial vector field. As they are orthogonal at each point, their Riemannian metric approximately vanishes.

# sample uniformly from a disc
n = 1200
disc_points =  np.random.uniform(-1, 1, (n, 2))
disc_points = disc_points[np.linalg.norm(disc_points, axis=1) <= 1]

# Create the diffusion geometry backbone
dg = DiffusionGeometry.from_point_cloud(disc_points, n_function_basis=200)

# Define a rotational vector field in pointwise basis
rot_pointwise = np.zeros_like(disc_points)
rot_pointwise[:, 0] = -disc_points[:, 1]
rot_pointwise[:, 1] = disc_points[:, 0]

# Construct VectorField object from pointwise basis representation
rot = VectorField.from_pointwise_basis(rot_pointwise, dg)

# Construct an outward radial vector field in pointwise basis
radial_pointwise = np.zeros_like(disc_points)
radial_pointwise[:, 0] = disc_points[:, 0]
radial_pointwise[:, 1] = disc_points[:, 1]
radial = VectorField.from_pointwise_basis(radial_pointwise, dg)

# Visualise the rotational vector field
scale = 0.12
line_width = 1

fig = make_subplots(
    rows=1,
    cols=3,
    specs=[[{"type": "xy"}]*3],
    horizontal_spacing=0.0,
    vertical_spacing=0.05,
    shared_xaxes=True,
    shared_yaxes=True,
    subplot_titles=("Rotational vector field: rot", "Radial vector field: radial", "Riemannian metric: g(rot, radial)")
)


fig.add_traces(plot_quiver_2d(dg.immersion_coords, rot.to_ambient(), scale=scale, line_width=line_width).data[0], rows=1, cols=1)
fig.add_traces(plot_quiver_2d(dg.immersion_coords, radial.to_ambient(), scale=scale, line_width=line_width).data[0], rows=1, cols=2)

# Compute the Riemannian metric between rot and radial
g_rot_radial = dg.g(rot, radial)

# This should approximate zero since the two vector fields are orthogonal
fig.add_traces(plot_scatter_2d(disc_points, color=g_rot_radial).data[0], rows=1, cols=3)

# We can quanitify this through the inner products:
print("Norm of rot:", rot.norm())
print("Norm of radial:", radial.norm())
print("Inner product of rot and radil:", dg.inner(rot, radial))

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

fig.update_layout(margin=dict(t=40))

Norm of rot: 0.6510319874450465
Norm of radial: 0.5665456686168123
Inner product of rot and radil: -0.0017776755967584615

Tangent Space Approximation

As one can see in the above plot, the ambient space representation of the vector field does not simply correspond to its initial pointwise description. More precisely, the vector fields are projections onto approximations of the tangent space. This tangent space inference is not explicit; it is a consequence of encoding the geometry via the Carré du champ and gets handled by the framework.

We can explicitly find a basis of the tangent space by forming the follow Gram matrix of ambient coordinate function: $$G(p)_{ij} := g_p(\nabla x_i(p), \nabla x_j(p))=\Gamma(x_i, x_j).$$ This matrix is symmetric. If the data is sampled from an $m$-dimensonal manifold, the $m$ eigenvectors with biggest eigenvalue serve as an approximation of a tangent basis. The other $d-m$ eigenvectors form a basis of the normal space.

We give the example of a circle and a torus.

# Tangent space approximation on the circle
# -----------------------------------------

n=50
xx = np.linspace(0, 2*np.pi, n)
circle_points = np.array([[np.cos(x), np.sin(x)] for x in xx])

dg_circle = DiffusionGeometry.from_point_cloud(circle_points, n_function_basis=n-1, n_coefficients=n-1)

# Access the Gram matrix
G = dg_circle.cache.gamma_coords

# Find the spectrum 
eigenvalues, eigenvectors = np.linalg.eigh(G)

tangent_vectors = eigenvectors[:, :, 1]
normal_vectors = eigenvectors[:, :, 0]



scale = 0.04
arrow_scale = 0.001

# Create subplots for tangent and normal vectors
fig = make_subplots(
    rows=1,
    cols=2,
    specs=[[{"type": "xy"},{"type": "xy"}]],
    horizontal_spacing=0.0,
    vertical_spacing=0.05,
    shared_xaxes=True,
    shared_yaxes=True,
    subplot_titles=("Tangent Vectors", "Normal Vectors")
)

# Draw tangent vectors and normal vectors
fig.add_traces(plot_quiver_2d(dg_circle.immersion_coords, tangent_vectors, scale=scale, arrow_scale=arrow_scale).data[0], rows=1, cols=1)
fig.add_traces(plot_quiver_2d(dg_circle.immersion_coords, -tangent_vectors, scale=scale, arrow_scale=arrow_scale).data[0], rows=1, cols=1)
fig.add_trace(plot_scatter_2d(dg_circle.immersion_coords, color='red').data[0])

fig.add_traces(plot_quiver_2d(dg_circle.immersion_coords, normal_vectors, scale=scale, arrow_scale=arrow_scale).data[0], rows=1, cols=2)
fig.add_traces(plot_quiver_2d(dg_circle.immersion_coords, -normal_vectors, scale=scale, arrow_scale=arrow_scale).data[0], rows=1, cols=2)
fig.add_traces(plot_scatter_2d(dg_circle.immersion_coords, color='red').data[0], rows=1, cols=2)

# Set equal aspect ratio
for r in range(1):
    for c in range(2):
        fig.update_xaxes(scaleanchor=f"y{c+1}", scaleratio=1, row=r+1, col=c+1)
        fig.update_yaxes(scaleanchor=f"x{c+1}", scaleratio=1, row=r+1, col=c+1)

clean_fig(fig)
fig.update_layout(margin=dict(t=40))

# Tangent space approximation on 3D spaces
# -----------------------------------------

# Code to generate surfaces 
# -----------------------------------------
def sample_sphere(n: int, radius: float = 1.0, seed=None) -> np.ndarray:
    if seed is not None:
        np.random.seed(seed)
    
    phi = np.random.uniform(0, 2*np.pi, n)       # azimuthal angle
    cos_theta = np.random.uniform(-1, 1, n)      # cosine of polar angle
    theta = np.arccos(cos_theta)                 # polar angle

    x = radius * np.sin(theta) * np.cos(phi)
    y = radius * np.sin(theta) * np.sin(phi)
    z = radius * cos_theta

    points = np.column_stack([x, y, z])
    return points

def sample_torus(n_points: int, R: float = 2.0, r: float = 0.7, seed=None) -> np.ndarray:
    if seed is not None:
        np.random.seed(seed)
    
    theta = np.random.uniform(0, 2*np.pi, n_points)  # tube angle
    phi = np.random.uniform(0, 2*np.pi, n_points)    # circle angle

    x = (R + r * np.cos(theta)) * np.cos(phi)
    y = (R + r * np.cos(theta)) * np.sin(phi)
    z = r * np.sin(theta)

    data = np.column_stack([x, y, z])
    return data, theta, phi

def sample_hyperboloid_one_sheet(n_u: int, n_v: int,
                                 a: float = 1.0, 
                                 b: float = 1.0,
                                 c: float = 1.0,
                                 u_max: float = 1.2,
                                 seed: int = None):

    if seed is not None:
        np.random.seed(seed)
    
    u_vals = np.random.uniform(-u_max, u_max, n_u * n_v)
    v_vals = np.random.uniform(0, 2*np.pi, n_u * n_v)

    x = a * np.cosh(u_vals) * np.cos(v_vals)
    y = b * np.cosh(u_vals) * np.sin(v_vals)
    z = c * np.sinh(u_vals)

    pts = np.column_stack([x, y, z])
    return pts, u_vals, v_vals


# ----------------------------------------------
# Extract tangent bundle and plot tangent planes
# ----------------------------------------------
def draw_tangent_planes(data, fig, column):
    dg = DiffusionGeometry.from_point_cloud(data, n_function_basis=200, n_coefficients=200)

    # Compute the tangent bundle
    _, eigenvectors = np.linalg.eigh(dg.cache.gamma_coords)

    # Since the surfaces are 2D, extract 2 eigenvectors with biggest eigenvalues
    v1 = eigenvectors[:, :, 1]
    v2 = eigenvectors[:, :, 2]
    bundle = np.stack((v1, v2), axis=2)

    # Plot tangent planes on the torus
    fig.add_traces(plot_tangent_planes_3d(bundle, data, angle=[1,1,1], zoom=1).data[0], rows=1, cols=column)

# Create subplot figure
fig = make_subplots(
    rows=1,
    cols=3,
    specs=[[{"type": "scene"}]*3],
    horizontal_spacing=0.0,
    vertical_spacing=0.05,
    shared_xaxes=True,
    shared_yaxes=True,
    subplot_titles=("Torus", "Sphere", "Hyperboloid")
)

# Draw tangent planes for each surface
torus_points, _, _ = sample_torus(2000)
draw_tangent_planes(torus_points, fig, column=1)

sphere_points = sample_sphere(2000)
draw_tangent_planes(sphere_points, fig, column=2)

hyperboloid_points, _, _ = sample_hyperboloid_one_sheet(40, 50)
draw_tangent_planes(hyperboloid_points, fig, column=3)


# Adjust figure layout
clean_fig(fig)
fig.update_layout(width=1000, height=400)
fig.update_layout(margin=dict(t=40))

fig.show()