Reproducing Theilman and Aimone 2025: a Neuromorphic Finite-Element Solver

PENG9570 Course Project Presentation

Michael Joseph Tarlton
Oslomet
June 2026

Introduction & Motivation

The Problem:

  • Modern PDE solver applications require ever larger levels of size and resolution
  • Modern PDE problems can easily require multiple-billions of individual points of resolution
  • In Numerical PDE solvers, e.g. the Finite Element Method (FEM), these points map to discretized individual elements in their overall integrative scheme
  • This is of course is extremely computationally expensive in conventional CPUs

The Neuromorphic Solution:

  • Analog computing may offer a solution
  • Spiking Neuromorphic Hardware (NMH) like Intel's Loihi 2 replicates neuronal dynamics directly in-silico
  • Continuous deterministic physical dynamics are modeled using discrete event-driven sparse spikes in a Spiking Neural Network (SNN) model
  • Extreme parallelization, high speed, and low power

The Challenge:

  • Implementing the FEM in an SNN is not a straightforward problem
  • Not only in mapping a FEM solver into an SNN structure
  • But also solving continuous valued problems in discretized spiking dynamics

Introducing NeuroFEM

Theilman and Aimone (2025) introduce an SNN implmentation of the FEM numerical solver [1]

NeuroFEM maps the discretized mesh of a FEM PDE solver directly into an SNN

  • Where each mesh node represents a discrete basis function in the FEM mesh, composed of a small local population of excitatory-inhibitory balanced integrate-and-fire neurons with recurrent connections
  • Mesh nodes are sparsely interconnected according to the mesh topology, e.g. tessellated across the surface on a disc

Introducing NeuroFEM

Theilman and Aimone (2025) introduce an SNN implmentation of the FEM numerical solver [1]

Driven dynamics instead of Training:

  • Does not learn via back-prop or iterative parameter search
  • The synaptic weights and are built once from the FEM matrix and the readout matrix , and then held fixed.
  • No learning rule modifies them. The per‑neuron state variables, integrated by forward Euler with evolves on every timestep
  • Collectively driving synaptic dynamics towards a solution
  • A forcing function correspondent to the function of the system is applied piecewise at each node, driving the local population of “fast” neurons to a local solution
  • While “slow” interconnected neurons integrate the activity from nearby nodes with the local population

Introducing NeuroFEM

Theilman and Aimone (2025) introduce an SNN implmentation of the FEM numerical solver [1]

An SNN implementation of a Proportional-Integral (PI) controller

  • The NeuroFEM builds on the balanced spiking network of Boerlin, Machens, and Denève (2013) [^2] which was Proportional only leaving a persistent steady-state offset
  • A proportional-integral (PI) controller drives a system's output toward a target by acting on the error between desired and actual values:

  • The proportional term responds to the present error.
  • The integral term accumulates past error. Its defining property: at steady state the control output can be nonzero even when , so the integral term can hold the system exactly on target and eliminate steady-state (offset) error that a proportional-only controller cannot.

Spiking Neural Network (LIF-PI) Dynamics

The Readout Variable ():
Nodal spiketrains are low-pass filtered through a balanced readout matrix to produce continuous state estimates:

Membrane Potential Dynamics ():

  • : slow synaptic integration (local estimate of )
  • : fast synaptic integration (within‑node coordination)
  • : local residual
    • Each neuron computes a local residual error in the “slow” synapse residual to zero using the running integral
  • : the running integral of
    • the PI controller's "I" term

Mathematical Methods: Poisson PDE

We solve the steady-state Poisson equation on a unit disk with Dirichlet boundary conditions:

The Weak Form (Galerkin Discretization):
Using piecewise-linear () triangular basis functions , we obtain the sparse linear system :

  • Symmetry & Positive Definiteness:
    The stiffness matrix is symmetric positive-definite, which is a mathematical requirement for stable recurrent SNN convergence.

Computational Implementation

  • Local Environment:
    • Evaluated end-to-end on conventional CPU using Python 3.12.
    • Notebooks: NeuroFEM_Annotated.ipynb and NeuroFEM_Sweep.ipynb.
    • Dependencies: NumPy, SciPy, Matplotlib, MeshPy/Triangle, tqdm.
  • Workflow Pipeline:
    1. Generate unstructured disk mesh using Shewchuk's Triangle generator.
    2. Assemble stiffness matrix and load vector via local element integration.
    3. Precondition the system and compile network weights (, , thresholds).
    4. Run forward-Euler integration of the LIF-PI equations for 50,000 timesteps ().

Baseline Results: Numerical Validation

Evaluated with interior mesh nodes, NPM neurons per node, and readout scaling . Time averaged over the final 10,000 steps.

Metric Description Value
Relative residual
NeuroFEM vs. Analytic Solution
NeuroFEM vs. SciPy Sparse Direct Solver
  • Analysis:

    • The difference between NeuroFEM and the direct solver is negligible
    • The NeuroFEM is a mathematically high-fidelity approximation of the physical continuous problem

Online Adaptability: Forcing Switch

  • The network is compiled once, but can resolve changing physical scenarios in real time.
  • We test this by switching the forcing bias at step 25,000:
    • (constant):
    • (quadratic):
Forcing switch result

Parameter Sweep: Resolution Scaling

We swept mesh counts up to 6503 and NPM values to reproduce the scaling claims in the original work.

Resolution scaling plot

Parameter Sweep: Analytic Error ()

  • Relative error decreases as the physical mesh is refined.
  • NeuroFEM (solid lines) tracks the exact conventional solver (black dashed) closely.
Analytic error plot

Parameter Sweep: System Residual ()

  • The residual per mesh node is independent of mesh scale
  • Improved directly by increasing NPM (local neural resolution) and decreasing .
System residual plot

Parameter Sweep: Solver Discrepancy ()

  • Absolute discrepancy between SNN output and conventional iterative solvers remains extremely small ().
Solver discrepancy plot

Limitations of Simulated Benchmarks

While our CPU simulation validates the mathematical framework of NeuroFEM, it omits key hardware-level advantages of real neuromorphic deployment:

  1. Energy & Power Constraints:
    • Loihi 2 leverages sparse, asynchronous event-driven pulses to consume microjoules per inference.
    • CPU simulations must calculate full floating-point equations for all nodes at every timestep, losing the physical energy scaling advantage.
  2. Fixed-Point Quantization:
    • Loihi 2 utilizes custom bit-width fixed-point integer values.
    • Our simulation assumes double-precision floating-point, hiding numerical precision bounds.

Conclusions & Next Steps

  • Conclusions:
    • Replicated the full NeuroFEM methodological chain.
    • Verified that PI controller successfully eliminates steady-state solution offsets.
    • Replicated baseline convergence, online forcing function transition, and parameter scaling sweeps.
  • Next Steps & Extensions:
    • Scale SNN models to 3D mechanical and elasticity problems.
    • Explore open-source neuromorphic simulators like the Lava Software Framework for closer hardware emulation.
    • Extend sparse mesh compilation paradigms to neural network architecture initialization.

References

[1] Theilman, B. H. & Aimone, J. B. (2025). Solving sparse finite element problems on neuromorphic hardware. Nature Machine Intelligence, 7(11), 1845-1857.
[2] Boerlin, M., Machens, C. K. & Denève, S. (2013). Predictive coding of dynamical variables in balanced spiking networks. PLoS Computational Biology, 9(11), e1003258.
[3] Davies, M. et al. (2021). Advancing neuromorphic computing with Loihi: A survey of results and outlook. Proceedings of the IEEE, 109(5), 911-934.
[4] Brenner, S. C. & Scott, L. R. (2008). The Mathematical Theory of Finite Element Methods, 3rd ed. Springer.

%% * Smith, J. D. et al. (2022). Neuromorphic scaling advantages for energy-efficient random walk computations. Nature Electronics, 5(2), 102-112. %%