21  Lab 21: Neuromorphic Computing: Building AI in the Brain’s Image

Open In Colab

21.1 Learning Objectives

  1. Understand core concepts from Chapter 21
  2. Implement key algorithms and techniques
  3. Apply methods to practical problems
  4. Analyze results and draw insights
  5. Connect to broader NeuroAI themes

21.2 Prerequisites

  • Reading: Chapter 21: Neuromorphic Computing
  • Libraries: NumPy, Matplotlib, relevant frameworks
  • Concepts: Event-driven processing, spiking neural networks

21.3 Setup

Run this cell to import the necessary libraries and configure the environment.

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(42)
plt.rcParams['figure.figsize'] = (12, 8)

21.4 Part 1: Event-Based Vision - Dynamic Vision Sensor Simulation

Core concepts and theory from Neuromorphic Computing. In this section, we implement a Dynamic Vision Sensor (DVS) simulation to understand how event-based cameras work.

21.4.1 Exercise 1: DVS Output Simulation

The DVS (Dynamic Vision Sensor) mimics the retina by generating events only when pixel brightness changes. Each event is a tuple (x, y, t, p) representing location, timestamp, and polarity (+1 for brightness increase, -1 for decrease).

Run this cell to implement the DVS simulation function.

def simulate_dvs_output(video_frames, threshold=0.1):
    """
    Simulate output of a Dynamic Vision Sensor from video frames.

    This function mimics how event-based cameras work: each pixel independently
    monitors brightness and generates an event when the log intensity changes
    by more than a threshold.

    Parameters
    ----------
    video_frames : ndarray
        Array of shape (T, H, W) containing T video frames
    threshold : float
        Log-intensity change threshold for triggering events

    Returns
    -------
    events : list of tuples
        Each event is (x, y, t, polarity) where:
        - x, y: pixel coordinates
        - t: timestamp (frame index)
        - polarity: +1 (brightness increase) or -1 (brightness decrease)
    """
    events = []
    prev_frame = video_frames[0]

    for t, frame in enumerate(video_frames[1:], 1):
        log_diff = np.log(frame + 1e-6) - np.log(prev_frame + 1e-6)
        on_events = np.where(log_diff > threshold)
        for y, x in zip(on_events[0], on_events[1]):
            events.append((x, y, t, 1))

        off_events = np.where(log_diff < -threshold)
        for y, x in zip(off_events[0], off_events[1]):
            events.append((x, y, t, -1))

        prev_frame = frame

    return events

21.4.2 Exercise 2: Generate Synthetic Video Frames

Create a sequence of synthetic video frames with a moving object to test the DVS simulation.

def generate_moving_bar_frames(num_frames=50, height=64, width=64, bar_width=5):
    """
    Generate video frames with a horizontal bar moving downward.

    Parameters
    ----------
    num_frames : int
        Number of frames to generate
    height, width : int
        Frame dimensions
    bar_width : int
        Width of the moving bar in pixels

    Returns
    -------
    frames : ndarray
        Array of shape (T, H, W) with pixel values in [0, 1]
    """
    frames = np.zeros((num_frames, height, width))

    for t in range(num_frames):
        # Bar position moves down over time
        bar_y = int((t / num_frames) * (height - bar_width))
        frames[t, bar_y:bar_y + bar_width, :] = 1.0

    return frames

# Generate test frames
video_frames = generate_moving_bar_frames(num_frames=50)
print(f"Generated video: {video_frames.shape}")

21.4.3 Exercise 3: Run DVS Simulation and Visualize Events

Apply the DVS simulation to the synthetic video and visualize the resulting events.

# Run DVS simulation
events = simulate_dvs_output(video_frames, threshold=0.1)
print(f"Generated {len(events)} events")

# Separate events by polarity
on_events = [(x, y, t) for x, y, t, p in events if p == 1]
off_events = [(x, y, t) for x, y, t, p in events if p == -1]

# Create visualization
fig, axes = plt.subplots(1, 3, figsize=(15, 5))

# Plot first frame
axes[0].imshow(video_frames[0], cmap='gray')
axes[0].set_title('Frame 0 (Start)')
axes[0].set_xlabel('x')
axes[0].set_ylabel('y')

# Plot last frame
axes[1].imshow(video_frames[-1], cmap='gray')
axes[1].set_title(f'Frame {len(video_frames)-1} (End)')
axes[1].set_xlabel('x')
axes[1].set_ylabel('y')

# Plot events as space-time diagram
if on_events:
    on_x, on_y, on_t = zip(*on_events)
    axes[2].scatter(on_x, on_t, c='green', s=1, alpha=0.5, label='ON events')
if off_events:
    off_x, off_y, off_t = zip(*off_events)
    axes[2].scatter(off_x, off_t, c='red', s=1, alpha=0.5, label='OFF events')

axes[2].set_xlabel('x position')
axes[2].set_ylabel('Time (frame)')
axes[2].set_title('Event Space-Time Diagram')
axes[2].legend()
axes[2].set_xlim(0, 64)
axes[2].set_ylim(0, 50)

plt.tight_layout()
plt.show()

21.4.4 Exercise 4: Analyze Data Reduction

Compare the amount of data generated by the DVS to conventional video frames.

# Calculate data reduction
total_pixels_per_frame = video_frames.shape[1] * video_frames.shape[2]  # H * W
total_frame_data = total_pixels_per_frame * len(video_frames)
event_data = len(events)

reduction_ratio = (1 - event_data / total_frame_data) * 100

print(f"\nData Reduction Analysis:")
print(f"  Frame data: {total_frame_data:,} pixel values")
print(f"  Event data: {event_data:,} events")
print(f"  Data reduction: {reduction_ratio:.1f}%")
print(f"  Compression ratio: {total_frame_data / event_data:.1f}x")

21.4.5 Exercise 5: Explore Threshold Sensitivity

Investigate how the threshold parameter affects event generation.

thresholds = [0.05, 0.1, 0.2, 0.5]
event_counts = []

for thresh in thresholds:
    evts = simulate_dvs_output(video_frames, threshold=thresh)
    event_counts.append(len(evts))
    print(f"Threshold {thresh}: {len(evts)} events")

# Plot threshold vs event count
fig, ax = plt.subplots(figsize=(8, 5))
ax.bar(range(len(thresholds)), event_counts, color='steelblue')
ax.set_xticks(range(len(thresholds)))
ax.set_xticklabels([str(t) for t in thresholds])
ax.set_xlabel('Threshold')
ax.set_ylabel('Number of Events')
ax.set_title('Effect of Threshold on Event Generation')
plt.tight_layout()
plt.show()

21.5 Part 2: Spiking Neural Network Fundamentals

21.5.1 Exercise 6: Leaky Integrate-and-Fire (LIF) Neuron

Implement a simple LIF neuron model to understand how spiking neurons process information.

def lif_neuron(I_input, dt=1.0, tau_m=20.0, v_threshold=1.0, v_reset=0.0, v_rest=0.0):
    """
    Leaky Integrate-and-Fire neuron model.

    Parameters
    ----------
    I_input : ndarray
        Input current over time
    dt : float
        Time step
    tau_m : float
        Membrane time constant
    v_threshold : float
        Spike threshold
    v_reset : float
        Reset potential after spike
    v_rest : float
        Resting potential

    Returns
    -------
    v_membrane : ndarray
        Membrane potential over time
    spikes : ndarray
        Binary spike train
    """
    n_steps = len(I_input)
    v_membrane = np.zeros(n_steps)
    spikes = np.zeros(n_steps)

    for t in range(1, n_steps):
        # Leaky integration
        dv = (-(v_membrane[t-1] - v_rest) + I_input[t-1]) / tau_m * dt
        v_membrane[t] = v_membrane[t-1] + dv

        # Spike generation
        if v_membrane[t] >= v_threshold:
            spikes[t] = 1
            v_membrane[t] = v_reset

    return v_membrane, spikes

# Test with step input
T = 200
I_input = np.zeros(T)
I_input[50:150] = 1.5  # Step current input

v_membrane, spikes = lif_neuron(I_input)

# Visualize
fig, axes = plt.subplots(3, 1, figsize=(12, 8), sharex=True)

axes[0].plot(I_input, 'b-', linewidth=2)
axes[0].set_ylabel('Input Current')
axes[0].set_title('LIF Neuron Response')

axes[1].plot(v_membrane, 'g-', linewidth=1.5)
axes[1].axhline(y=1.0, color='r', linestyle='--', label='Threshold')
axes[1].set_ylabel('Membrane Potential')
axes[1].legend()

axes[2].stem(np.where(spikes)[0], spikes[spikes == 1], basefmt=' ')
axes[2].set_ylabel('Spikes')
axes[2].set_xlabel('Time (ms)')

plt.tight_layout()
plt.show()

print(f"Total spikes: {int(spikes.sum())}")
print(f"Spike times: {np.where(spikes)[0]}")

21.5.2 Exercise 7: Spiking Network Layer

Create a small network of spiking neurons to observe population dynamics.

def spiking_layer(input_currents, weights, n_neurons=10, dt=1.0, tau_m=20.0):
    """
    A layer of spiking neurons with random recurrent connections.

    Parameters
    ----------
    input_currents : ndarray
        Input current over time (T,)
    weights : ndarray
        Weight matrix (n_neurons, n_neurons) for recurrent connections
    n_neurons : int
        Number of neurons in the layer
    dt, tau_m : float
        Time step and membrane time constant

    Returns
    -------
    spikes : ndarray
        Spike raster (T, n_neurons)
    """
    T = len(input_currents)
    v_membrane = np.zeros(n_neurons)
    spikes = np.zeros((T, n_neurons))

    v_threshold = 1.0
    v_reset = 0.0

    for t in range(T):
        # Input plus recurrent input from previous spikes
        if t > 0:
            recurrent_input = weights @ spikes[t-1]
        else:
            recurrent_input = np.zeros(n_neurons)

        I_total = input_currents[t] + recurrent_input

        # Update membrane potential
        dv = (-v_membrane + I_total) / tau_m * dt
        v_membrane = v_membrane + dv

        # Generate spikes
        fired = v_membrane >= v_threshold
        spikes[t] = fired.astype(float)
        v_membrane[fired] = v_reset

    return spikes

# Create network
n_neurons = 10
np.random.seed(42)
weights = np.random.randn(n_neurons, n_neurons) * 0.3
np.fill_diagonal(weights, 0)  # No self-connections

# Generate input
T = 300
input_currents = np.zeros(T)
input_currents[50:250] = 0.8 + 0.3 * np.sin(np.linspace(0, 4*np.pi, 200))

# Run simulation
spikes = spiking_layer(input_currents, weights, n_neurons)

# Visualize raster plot
fig, ax = plt.subplots(figsize=(12, 6))

for neuron_idx in range(n_neurons):
    spike_times = np.where(spikes[:, neuron_idx])[0]
    ax.scatter(spike_times, [neuron_idx] * len(spike_times),
               marker='|', s=100, c='black')

ax.set_xlabel('Time (ms)')
ax.set_ylabel('Neuron Index')
ax.set_title('Spiking Network Raster Plot')
ax.set_ylim(-0.5, n_neurons - 0.5)

plt.tight_layout()
plt.show()

# Calculate sparsity
total_possible_spikes = T * n_neurons
actual_spikes = spikes.sum()
sparsity = 1 - (actual_spikes / total_possible_spikes)

print(f"\nSparsity Analysis:")
print(f"  Total possible spikes: {total_possible_spikes}")
print(f"  Actual spikes: {int(actual_spikes)}")
print(f"  Sparsity: {sparsity * 100:.1f}%")

21.6 Part 3: Memristor Crossbar Array Simulation

21.6.1 Exercise 8: In-Memory Computing

Simulate a memristor crossbar array for matrix-vector multiplication, demonstrating the principle of in-memory computing.

def memristor_crossbar_mvm(input_voltages, conductance_matrix):
    """
    Simulate memristor crossbar array performing matrix-vector multiplication.

    In a crossbar array, input voltages are applied to rows, and each memristor
    at the intersection encodes a weight as its conductance. Output currents
    at columns are computed via Ohm's law (I = V * G).

    Parameters
    ----------
    input_voltages : ndarray
        Input vector (n_inputs,)
    conductance_matrix : ndarray
        Memristor conductance values (n_outputs, n_inputs)

    Returns
    -------
    output_currents : ndarray
        Output currents from columns (n_outputs,)
    """
    # Memristor crossbar performs: output = conductance @ input
    # This is naturally parallel in hardware
    output_currents = conductance_matrix @ input_voltages
    return output_currents

# Create a 4x4 crossbar array
np.random.seed(42)
n_inputs = 4
n_outputs = 4

# Memristor conductances (normalized for simulation)
conductance_matrix = np.random.rand(n_outputs, n_inputs) * 0.5 + 0.1

# Input voltages
input_voltages = np.array([0.5, 0.8, 0.2, 0.9])

# Compute via crossbar (parallel in hardware)
output_currents = memristor_crossbar_mvm(input_voltages, conductance_matrix)

# Verify against standard matrix multiplication
expected_output = conductance_matrix @ input_voltages

print("Memristor Crossbar Array Simulation")
print("=" * 40)
print(f"\nInput voltages: {input_voltages}")
print(f"\nConductance matrix:\n{conductance_matrix}")
print(f"\nOutput currents (crossbar): {output_currents}")
print(f"Expected (standard MVM):    {expected_output}")
print(f"\nVerification: {np.allclose(output_currents, expected_output)}")

# Compare computational cost
standard_ops = n_inputs * n_outputs  # Multiply-accumulate operations
crossbar_ops = 1  # Single parallel step in hardware

print(f"\nComputational Cost Comparison:")
print(f"  Standard MVM: {standard_ops} sequential multiply-accumulate ops")
print(f"  Crossbar: {crossbar_ops} parallel step (all {standard_ops} MACs simultaneously)")
print(f"  Hardware speedup: {standard_ops}x")

21.7 Part 4: Energy-Delay Product Analysis

21.7.1 Exercise 9: Efficiency Metrics

Analyze the Energy-Delay Product (EDP) to understand neuromorphic efficiency advantages.

# Hypothetical benchmark data
platforms = {
    'GPU (V100)': {'energy_J': 5000, 'latency_ms': 50},
    'CPU (Xeon)': {'energy_J': 8000, 'latency_ms': 100},
    'Loihi 2': {'energy_J': 5, 'latency_ms': 20},
    'TrueNorth': {'energy_J': 2, 'latency_ms': 10}
}

# Calculate EDP for each platform
print("Energy-Delay Product Analysis")
print("=" * 50)
print(f"{'Platform':<15} {'Energy (J)':<12} {'Latency (ms)':<14} {'EDP':<15}")
print("-" * 50)

edp_values = {}
for name, specs in platforms.items():
    edp = specs['energy_J'] * specs['latency_ms']
    edp_values[name] = edp
    print(f"{name:<15} {specs['energy_J']:<12} {specs['latency_ms']:<14} {edp:<15,}")

print("-" * 50)

# Calculate improvement ratios
gpu_edp = edp_values['GPU (V100)']
print(f"\nImprovement vs GPU:")
for name, edp in edp_values.items():
    improvement = gpu_edp / edp
    print(f"  {name}: {improvement:.0f}x better EDP")

# Visualization
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# Energy vs Latency scatter plot
names = list(platforms.keys())
energies = [platforms[n]['energy_J'] for n in names]
latencies = [platforms[n]['latency_ms'] for n in names]
colors = ['red', 'orange', 'blue', 'green']

for i, name in enumerate(names):
    axes[0].scatter(latencies[i], energies[i], c=colors[i], s=200, label=name, marker='o')

axes[0].set_xlabel('Latency (ms)')
axes[0].set_ylabel('Energy (J)')
axes[0].set_title('Energy vs Latency Trade-off')
axes[0].legend()
axes[0].set_yscale('log')
axes[0].grid(True, alpha=0.3)

# EDP bar chart
bars = axes[1].bar(range(len(names)), [edp_values[n] for n in names], color=colors)
axes[1].set_xticks(range(len(names)))
axes[1].set_xticklabels(names, rotation=15)
axes[1].set_ylabel('Energy-Delay Product (J*ms)')
axes[1].set_title('Energy-Delay Product Comparison')
axes[1].set_yscale('log')
axes[1].grid(True, alpha=0.3, axis='y')

plt.tight_layout()
plt.show()

21.8 Exercises

21.8.1 Exercise 10: Implement Your Own DVS Test

Create your own synthetic video scenario and analyze the DVS output.

def generate_expanding_circle_frames(num_frames=60, height=64, width=64):
    """
    Generate video frames with an expanding circle.

    TODO: Implement this function to create frames where a circle
    centered at (height/2, width/2) expands over time.

    Hint: Use np.ogrid to create coordinate arrays and check if
    distance from center is less than the current radius.
    """
    # Your implementation here
    frames = np.zeros((num_frames, height, width))

    # Create coordinate grids
    y, x = np.ogrid[:height, :width]
    center_y, center_x = height // 2, width // 2

    for t in range(num_frames):
        # Radius grows from 5 to 25 pixels
        radius = 5 + int((t / num_frames) * 20)
        distance = np.sqrt((x - center_x)**2 + (y - center_y)**2)
        frames[t] = (distance < radius).astype(float)

    return frames

# Test your implementation
circle_frames = generate_expanding_circle_frames()
events = simulate_dvs_output(circle_frames, threshold=0.1)
print(f"Expanding circle generated {len(events)} events")

21.8.2 Exercise 11

Explore parameter variations in the LIF neuron model.

21.8.3 Exercise 12

Apply to real or simulated data.

21.8.4 Exercise 13

Compare different approaches for encoding information in spikes.

21.8.5 Exercise 14

Discuss NeuroAI connections between biological and artificial spiking systems.


21.9 Challenge Problems

21.9.1 Challenge 1

Implement a complete SNN for a classification task using surrogate gradient training.

21.9.2 Challenge 2

Reproduce research findings from a neuromorphic computing paper (e.g., event-based object detection).

21.9.3 Challenge 3

Design a novel application combining event cameras with spiking neural networks.


Open In Colab

21.10 Discussion Questions

21.10.1 Question 1

How does the event-based paradigm of DVS cameras relate to the way biological retinas process visual information?

21.10.2 Question 2

What are the trade-offs between rate coding (more spikes = stronger signal) and temporal coding (precise spike timing matters) in SNNs?

21.10.3 Question 3

How might on-chip learning with STDP enable truly autonomous systems that can adapt to their environment?


21.11 Summary

Explored Event-driven processing and neuromorphic computing through theory and practice.

Key Takeaways: - Dynamic Vision Sensors generate sparse, event-based data with dramatic data reduction - Spiking Neural Networks achieve efficiency through sparse, event-driven computation - Memristor crossbar arrays enable in-memory computing, eliminating the von Neumann bottleneck - Neuromorphic systems offer orders-of-magnitude improvements in Energy-Delay Product for sparse workloads - Core concepts mastered through hands-on implementation - Practical implementation skills developed - NeuroAI connections understood between biological and artificial systems