discrete-heat-kernels-simplicial

name	discrete-heat-kernels-simplicial
description	Discrete Heat Kernels on Simplicial Complexes and Its Application to Functional Brain Networks. Unified framework for heat kernel smoothing on simplicial complexes extending classical signal processing methods to higher-order network structures.

Discrete Heat Kernels on Simplicial Complexes

Unified framework for heat kernel smoothing on simplicial complexes, extending classical signal processing methods to higher-order network structures for functional brain network analysis.

Metadata

• Source: arXiv:2509.16908v2
• Authors: Sixtus Dakurah
• Published: 2025-09-21
• Category: cs.SI, q-bio.NC, math.CO

Core Methodology

Key Innovation

Conventional graph-theoretic analyses capture exclusively pairwise interactions, omitting intricate higher-order relationships. This methodology introduces a unified framework for heat kernel smoothing on simplicial complexes, extending classical signal processing to higher-order structures.

Theoretical Foundation

Heat Equation on Graphs (Review)

Standard Graph Heat Equation:

∂u/∂t = -L u

where L = D - A (Laplacian)
      u(t) = vertex signal at time t

Solution - Heat Kernel:

u(t) = e^(-tL) u(0)

Heat kernel H_t = e^(-tL) = Σ_k e^(-tλ_k) φ_k φ_k^T

where λ_k, φ_k are eigenvalues/vectors of L

Properties:

• Smoothing increases with time t
• Eigenvalues λ_k control frequency
• High λ_k (high frequency) decay faster

Extension to Simplicial Complexes

Simplicial Complex Structure:

• k-simplices: {v_0, ..., v_k} where (k+1) nodes all interact
• Boundary operator ∂_k: maps k-simplices to (k-1)-simplices
• Coboundary operator δ_k: adjoint of ∂_k

Hodge Laplacian:

L_k = ∂_{k+1} ∘ δ_k + δ_{k-1} ∘ ∂_k

= L_k^{up} + L_k^{down}

where L_k^{up} = ∂_{k+1} δ_k      (upper Laplacian)
      L_k^{down} = δ_{k-1} ∂_k    (lower Laplacian)

Interpretation:

• L_0: Standard graph Laplacian (on vertices)
• L_1: Edge Laplacian - captures flow/circulation
• L_2: Face Laplacian - captures voids/cavities

Heat Equation on k-simplices

∂u_k/∂t = -L_k u_k

Solution: u_k(t) = e^(-t L_k) u_k(0)

Key Insight: Heat diffusion occurs simultaneously on all dimensions, capturing different aspects of network dynamics.

Heat Kernel Applications

1. Signal Smoothing

Standard (Vertex-level):

u_smoothed = H_t^(0) * u_0

Higher-order (Edge-level):

flow_smoothed = H_t^(1) * f_0

where f_0 is initial edge flow

Application to fMRI:

• Vertex signals: Regional BOLD activity
• Edge signals: Information flow between regions
• Face signals: Multi-region coordination

2. Distance Metrics

Diffusion Distance:

d_t(i, j)^2 = ||H_t(δ_i - δ_j)||^2

= Σ_k e^(-2tλ_k) (φ_k(i) - φ_k(j))^2

Interpretation: Measures connectivity strength considering all paths, not just shortest.

3. Clustering and Segmentation

Spectral Clustering on SC:

1. Compute L_k for target dimension k
2. Eigen-decomposition: L_k = Φ Λ Φ^T
3. Use first m eigenvectors for k-means
4. Clusters represent cohesive higher-order structures

Functional Brain Network Pipeline

Preprocessed fMRI
    ↓
1. Time Series Extraction
    - ROI-based parcellation
    - BOLD signal extraction
    ↓
2. Simplicial Complex Construction
    - Compute correlations
    - Threshold to create edges
    - Find cliques → simplices
    ↓
3. Heat Kernel Computation
    - Build Hodge Laplacian L_k
    - Compute eigen-decomposition
    - Construct H_t = e^(-tL_k)
    ↓
4. Signal Processing
    - Apply heat kernel smoothing
    - Extract features
    - Classification/Analysis

Implementation Guide

Prerequisites

# Required libraries
pip install numpy scipy
pip install scikit-learn
pip install networkx
pip install gudhi          # For persistent homology/simplicial complexes
pip install nilearn        # fMRI processing

Step-by-Step Implementation

Step 1: Simplicial Complex Construction from fMRI

import numpy as np
from itertools import combinations
import networkx as nx

def build_simplicial_complex(time_series, correlation_threshold=0.5, max_dim=3):
    """
    Build simplicial complex from fMRI time series.
    
    Args:
        time_series: [n_timepoints, n_regions] array
        correlation_threshold: Minimum correlation for edge
        max_dim: Maximum simplex dimension
    
    Returns:
        simplices: List of (dimension, vertices) tuples
        L_k: Dictionary of Hodge Laplacians
    """
    n_regions = time_series.shape[1]
    
    # Compute correlation matrix
    corr_matrix = np.corrcoef(time_series.T)
    
    # Build graph from thresholded correlations
    G = nx.Graph()
    G.add_nodes_from(range(n_regions))
    
    edges = []
    for i in range(n_regions):
        for j in range(i+1, n_regions):
            if abs(corr_matrix[i, j]) > correlation_threshold:
                G.add_edge(i, j, weight=corr_matrix[i, j])
                edges.append((i, j))
    
    # Find cliques and build simplices
    simplices = {0: [(i,) for i in range(n_regions)]}
    
    for k in range(1, max_dim + 1):
        simplices[k] = []
        # Find all (k+1)-cliques
        for clique in nx.find_cliques(G):
            if len(clique) >= k + :
                 simplex  combinations(clique, k + ):
                    simplices[k].append(((simplex)))
        
        simplices[k] = ((simplices[k]))
    
     simplices, G, corr_matrix


time_series = np.random.randn(, )  
simplices, graph, corr = build_simplicial_complex(time_series, max_dim=)

()
()
()

Step 2: Boundary Operators and Hodge Laplacian

def construct_boundary_operators(simplices):
    """
    Construct boundary operators ∂_k for k = 1, 2, ..., max_dim.
    
    Args:
        simplices: Dictionary {k: list of k-simplices}
    
    Returns:
        boundary_ops: Dictionary {k: ∂_k as sparse matrix}
    """
    from scipy.sparse import csr_matrix
    
    boundary_ops = {}
    max_dim = max(simplices.keys())
    
    for k in range(1, max_dim + 1):
        k_simplices = simplices[k]
        k_minus_1_simplices = simplices[k-1]
        
        # Create index mapping for (k-1)-simplices
        idx_map = {simplex: i for i, simplex in enumerate(k_minus_1_simplices)}
        
        rows = []
        cols = []
        data = []
        
        for j, k_simplex in enumerate(k_simplices):
            # Faces of k-simplex (oriented boundary)
            for i, vertex in enumerate(k_simplex):
                face = k_simplex[:i] + k_simplex[i+1:]
                orientation = (-1) ** i
                
                if face in idx_map:
                    rows.append(idx_map[face])
                    cols.append(j)
                    data.append(orientation)
        
        boundary_ops[k] = csr_matrix(
            (data, (rows, cols)),
            shape=(len(k_minus_1_simplices), len(k_simplices))
        )
    
    return boundary_ops


def compute_hodge_laplacian(boundary_ops, k):
    """
    Compute Hodge Laplacian L_k.
    
    L_k = ∂_{k+1} ∂_{k+1}^T + ∂_k^T ∂_k
        = L_k^{up} + L_k^{down}
    """
     k == :
        
        B1 = boundary_ops[]
        L_0 = B1 @ B1.T
         L_0
    
    B_k = boundary_ops[k]
    L_down = B_k.T @ B_k
    
     k +   boundary_ops:
        B_k_plus_1 = boundary_ops[k + ]
        L_up = B_k_plus_1 @ B_k_plus_1.T
    :
        L_up = 
    
    L_k = L_up + L_down
     L_k



boundary_ops = construct_boundary_operators(simplices)
L_0 = compute_hodge_laplacian(boundary_ops, )  
L_1 = compute_hodge_laplacian(boundary_ops, )  
L_2 = compute_hodge_laplacian(boundary_ops, )  

()
()
()

Step 3: Heat Kernel Computation

def compute_heat_kernel(L_k, time_points, method='eigen'):
    """
    Compute heat kernel H_t = e^(-t * L_k).
    
    Args:
        L_k: Hodge Laplacian (sparse or dense matrix)
        time_points: List of time values t
        method: 'eigen' (eigendecomposition) or 'taylor' (series expansion)
    
    Returns:
        heat_kernels: Dictionary {t: H_t}
    """
    from scipy.linalg import expm
    from scipy.sparse.linalg import eigsh
    
    n = L_k.shape[0]
    heat_kernels = {}
    
    if method == 'eigen':
        # Eigendecomposition (for smaller matrices)
        # For sparse matrices, use eigsh for partial decomposition
        if n < 1000:
            eigenvalues, eigenvectors = np.linalg.eigh(L_k.toarray() if hasattr(L_k, 'toarray') else L_k)
        else:
            # Partial eigendecomposition for large matrices
            k = min(n - 2, 100)  # Top k eigenvalues
            eigenvalues, eigenvectors = eigsh(L_k, k=k, which='SM')
        
        for t in time_points:
            # H_t = Φ * diag(e^(-t*λ)) * Φ^T
            H_t = eigenvectors @ np.diag(np.exp(-t * eigenvalues)) @ eigenvectors.T
            heat_kernels[t] = H_t
    
    elif method == 'taylor':
        # Matrix exponential via scipy
        L_dense = L_k.toarray() if hasattr(L_k, 'toarray') else L_k
        for t in time_points:
            H_t = expm(-t * L_dense)
            heat_kernels[t] = H_t
    
     heat_kernels



time_points = [, , , , ]
H_0 = compute_heat_kernel(L_0, time_points)  
H_1 = compute_heat_kernel(L_1, time_points)  

()
()

Step 4: Signal Smoothing and Feature Extraction

class SimplicialSignalProcessor:
    """Process signals on simplicial complexes using heat kernels."""
    
    def __init__(self, simplices, boundary_ops):
        self.simplices = simplices
        self.boundary_ops = boundary_ops
        self.heat_kernels = {}
    
    def compute_all_heat_kernels(self, time_points):
        """Pre-compute heat kernels for all dimensions."""
        for k in self.simplices.keys():
            if k < max(self.simplices.keys()):
                L_k = compute_hodge_laplacian(self.boundary_ops, k)
                self.heat_kernels[k] = compute_heat_kernel(L_k, time_points)
    
    def smooth_signal(self, signal, k, t):
        """
        Apply heat kernel smoothing to k-dimensional signal.
        
        Args:
            signal: Initial signal on k-simplices [n_k]
            k: Simplex dimension
            t: Diffusion time
        
        Returns:
            smoothed: Smoothed signal
        """
        H_t = self.heat_kernels[k][t]
        smoothed = H_t @ signal
        return smoothed
    
    def extract_features(self, time_series, time_points):
        """
        Extract heat kernel features for classification.
        
        Args:
            time_series: Brain activity time series
            time_points: List of diffusion times
        
        Returns:
            features: Feature vector
        """
        features = {}
        
        # Vertex-level features (standard BOLD signals)
        v_signal = time_series.mean(axis=)  
         t  time_points:
            v_smooth = .smooth_signal(v_signal, , t)
            features[] = np.mean(v_smooth)
            features[] = np.std(v_smooth)
        
        
           .heat_kernels:
            
            edge_signal = np.array([
                corr[e]  e  .simplices[]
            ])
             t  time_points:
                e_smooth = .smooth_signal(edge_signal, , t)
                features[] = np.mean(e_smooth)
                features[] = np.(e_smooth**)
        
         features



processor = SimplicialSignalProcessor(simplices, boundary_ops)
processor.compute_all_heat_kernels([, , , ])


features = processor.extract_features(time_series, [, , , ])
()

Step 5: Complete fMRI Analysis Pipeline

class fMRIHeatKernelAnalysis:
    """Complete pipeline for fMRI analysis using heat kernels on SC."""
    
    def __init__(self, max_dim=2, time_points=None):
        self.max_dim = max_dim
        self.time_points = time_points or [0.1, 0.5, 1.0, 2.0, 5.0]
    
    def process_subject(self, fmri_data, atlas):
        """
        Process single subject fMRI data.
        
        Args:
            fmri_data: 4D fMRI volume
            atlas: Parcellation atlas
        
        Returns:
            features: Dictionary of features
        """
        # Extract time series
        from nilearn.maskers import NiftiLabelsMasker
        masker = NiftiLabelsMasker(labels_img=atlas, standardize=True)
        time_series = masker.fit_transform(fmri_data)
        
        # Build simplicial complex
        simplices, graph, corr = build_simplicial_complex(
            time_series, 
            max_dim=self.max_dim
        )
        
        # Compute boundary operators and heat kernels
        boundary_ops = construct_boundary_operators(simplices)
        processor = SimplicialSignalProcessor(simplices, boundary_ops)
        processor.compute_all_heat_kernels(self.time_points)
        
        # Extract features
        features = processor.extract_features(time_series, self.time_points)
        
        return features
    
    def analyze_group(self, subject_list, labels=None):
        
        all_features = []
        
         fmri_file, atlas  subject_list:
            fmri_data = nibabel.load(fmri_file)
            features = .process_subject(fmri_data, atlas)
            all_features.append(features)
        
        
         pandas  pd
        df = pd.DataFrame(all_features)
        
         df

Applications

• Functional connectivity smoothing with higher-order context
• Network distance metrics for brain connectivity analysis
• Multi-scale clustering of brain regions
• Disease classification (ADHD, ASD, Alzheimer's)
• Dynamic network analysis with time-varying heat kernels

Advantages Over Standard Graph Methods

Aspect	Graph Laplacian	Simplicial Heat Kernel
Structure	Vertices, edges	Vertices, edges, faces, ...
Laplacian	L_0 only	L_0, L_1, L_2, ...
Signal	Vertex signals	Signals on all dimensions
Smoothing	Local averaging	Multi-scale diffusion
Higher-order	Not captured	Explicitly modeled

Pitfalls

• Computational cost: Hodge Laplacian computation scales with simplex count
• Simplex explosion: Large cliques create many higher-order simplices
• Threshold sensitivity: Correlation threshold affects simplex structure
• Interpretation: Higher-order features need domain expertise
• Eigenvalue computation: Large matrices require approximate methods

Related Skills

• combinatorial-complex-brain-fmri
• higher-order-brain-networks
• dcho-higher-order-brain-connectivity
• brain-higher-order-structures

References

@article{dakurah2025heatkernels,
  title={Discrete Heat Kernels on Simplicial Complexes and Its Application to Functional Brain Networks},
  author={Dakurah, Sixtus},
  journal={arXiv preprint arXiv:2509.16908},
  year={2025}
}