Micrologics Logo
Back to Blog
AIPublished on July 20, 2026

The Architecture of Probabilistic Computing: Simulating p-Bits in Software for Combinatorial Optimization

Discover how probabilistic computing leverages thermal noise to solve NP-hard optimization problems. Learn the mathematical foundations of p-bits and how to build a high-performance vectorised simulator in Python.

Beyond Determinism: The Rise of the Probabilistic Computer

For decades, the trajectory of computer architecture has been defined by a relentless pursuit of determinism. We design transistors to switch cleanly between 0 and 1, build error-correcting memory to suppress thermal fluctuations, and write compilers that guarantee predictable execution paths. Yet, as we push against the physical limits of silicon and confront NP-hard optimization problems, this deterministic paradigm is becoming a bottleneck.

Enter probabilistic computing. Unlike classical computers that treat noise as an enemy to be eliminated, probabilistic computers leverage thermal noise as a computational resource. By replacing standard bits with probabilistic bits (p-bits), we can solve complex combinatorial optimization tasks, simulate quantum-like systems at room temperature, and accelerate machine learning algorithms at a fraction of the energy cost.

In this article, we will dissect the architecture of probabilistic computing, explore the mathematics behind p-bits, and build a software-defined simulator from scratch to solve a classic optimization problem.

The Anatomy of a Probabilistic Bit (p-Bit)

To understand a p-bit, it is helpful to compare it to both a classical bit and a quantum bit (qubit):

  • Classical Bit: Holds a static state of either 0 or 1.
  • Qubit: Exists in a coherent superposition of |0> and |1>, collapsing to a deterministic state only upon measurement. Qubits require extreme cryogenic temperatures to maintain coherence.
  • p-Bit: A fluctuating system that rapidly transitions between 0 and 1 over time. The probability of finding the p-bit in state 1 or 0 is governed by an input bias signal.

Mathematically, the state of a p-bit $s_i \in {-1, 1}$ at time $t$ is determined by an input $I_i(t)$ through a non-linear activation function, typically a hyperbolic tangent, modulated by thermal noise. The probability $P(s_i = 1)$ is expressed as:

$$P(s_i(t) = 1) = \frac{1}{2} [ 1 + \tanh(\beta \cdot I_i(t)) ]$$

Where:

  • $\beta$ represents the inverse temperature (or noise control parameter).
  • $I_i(t)$ is the input bias, which is a weighted sum of the states of other connected p-bits plus an external bias $h_i$.

If the input $I_i$ is highly positive, the p-bit is almost certainly in state 1. If it is highly negative, it is in state -1. If $I_i$ is zero, the p-bit fluctuates randomly between -1 and 1 with equal probability.

The Synaptic Network: Solving Combinatorial Optimization

To solve problems with p-bits, we connect them in a network where the coupling weights represent the constraints of our problem. This is mathematically isomorphic to an Ising Spin Glass model or a Quadratic Unconstrained Binary Optimization (QUBO) problem.

The total energy of the system is defined by the Hamiltonian:

$$E(s) = -\sum_{i < j} J_{ij} s_i s_j - \sum_i h_i s_i$$

Where $J_{ij}$ is the coupling strength between p-bit $i$ and p-bit $j$, and $h_i$ is the local bias of p-bit $i$. The system naturally gravitates toward states of lower energy. By letting the p-bits fluctuate and slowly cooling the system (increasing $\beta$), the network settles into the global minimum of the energy landscape, which corresponds to the optimal solution of our problem.

Building a Software-Based p-Bit Simulator

While dedicated hardware (such as magnetoresistive random-access memory or MRAM-based p-bits) is the ultimate goal, we can simulate a probabilistic computer in software to understand its dynamics.

Let us implement a vectorized p-bit simulator in Python using NumPy to solve a Max-Cut problem. The Max-Cut problem asks us to divide the vertices of a graph into two sets such that the number of edges between the two sets is maximized.

import numpy as np

class ProbabilisticComputer:
    def __init__(self, num_pbits, coupling_matrix, local_biases=None):
        self.N = num_pbits
        self.J = coupling_matrix  # Symmetric matrix of size (N, N)
        self.h = local_biases if local_biases is not None else np.zeros(self.N)
        # Initialize p-bits randomly to -1 or 1
        self.states = np.random.choice([-1, 1], size=self.N).astype(np.float64)

    def step(self, beta):
        """
        Performs a single update cycle across all p-bits.
        For a true parallel simulation, we update asynchronously or use a color-coding scheme.
        """
        # Calculate the input to each p-bit: I = J * s + h
        inputs = np.dot(self.J, self.states) + self.h

        # Calculate the probability of being in state +1
        probabilities = 0.5 * (1.0 + np.tanh(beta * inputs))

        # Generate random noise
        random_draws = np.random.rand(self.N)

        # Update states based on the probability threshold
        self.states = np.where(random_draws < probabilities, 1.0, -1.0)

    def get_energy(self):
        """Calculates the current Ising energy of the system."""
        interaction_energy = -0.5 * np.dot(self.states, np.dot(self.J, self.states))
        bias_energy = -np.dot(self.h, self.states)
        return interaction_energy + bias_energy

Simulating the Annealing Process

To find the optimal state, we must implement an annealing schedule. We start with a low $\beta$ (high temperature/high noise) to allow the system to explore the state space freely, and gradually increase $\beta$ (low temperature/low noise) to freeze the system into its lowest-energy state.

def solve_max_cut(graph_adjacency, steps=1000, beta_start=0.1, beta_end=3.0):
    # For Max-Cut, J_ij = -W_ij (negative of edge weights)
    J = -graph_adjacency
    num_nodes = graph_adjacency.shape[0]

    computer = ProbabilisticComputer(num_pbits=num_nodes, coupling_matrix=J)

    best_states = np.copy(computer.states)
    best_energy = computer.get_energy()

    # Linear annealing schedule
    betas = np.linspace(beta_start, beta_end, steps)

    for step, beta in enumerate(betas):
        computer.step(beta)
        current_energy = computer.get_energy()

        if current_energy < best_energy:
            best_energy = current_energy
            best_states = np.copy(computer.states)

    return best_states, best_energy

Optimizing the Simulator for Scale

Running sequential updates on a CPU is highly inefficient for large-scale p-bit networks. To simulate thousands or millions of p-bits, we must address two primary bottlenecks: parallel updates and fast random number generation.

1. Asynchronous vs. Synchronous Updates

If we update all p-bits simultaneously (synchronously), the system can fall into chaotic oscillations. For example, if two strongly coupled p-bits both see that the other is $+1$, they might both switch to $-1$ in the next step, and then back to $+1$, oscillating indefinitely.

To prevent this, we must use:

  • Asynchronous Updates: Update one random p-bit at a time. This is mathematically rigorous but slow to simulate.
  • Graph Coloring: Partition the graph of p-bits into independent sets (colors) where no two connected nodes share the same color. We can then update all nodes of a single color in parallel without race conditions or oscillations.

2. High-Throughput PRNGs

The bottleneck in simulating probabilistic systems is almost always the Pseudo-Random Number Generator (PRNG). Standard generators like Mersenne Twister are too slow for high-throughput execution.

For high-performance GPU-accelerated simulations (using PyTorch or CUDA), we use lightweight generators like Xorshift or PCG (Permuted Congruential Generator) implemented directly in the execution kernels.

Here is how you can write a highly vectorized PyTorch implementation that leverages the GPU for massive parallelism:

import torch

def torch_pbit_update(states, J, h, beta):
    # states: (batch_size, N)
    # J: (N, N)
    # h: (N,)
    inputs = torch.matmul(states, J.t()) + h
    probabilities = 0.5 * (1.0 + torch.tanh(beta * inputs))
    random_draws = torch.rand_like(states)
    return torch.where(random_draws < probabilities, 1.0, -1.0)

Using batching, we can run thousands of independent annealing paths in parallel, drastically increasing the probability of finding the global minimum in highly complex energy landscapes.

The Future: Physical Probabilistic Hardware

While software simulation is invaluable for algorithm design, the true potential of probabilistic computing lies in dedicated hardware. Silicon-based p-bits can be built using unstable magnetic tunnel junctions (MTJs)—the same technology used in MRAM. By intentionally designing MTJs to be thermally unstable, they switch states at gigahertz frequencies, providing a natural source of high-speed, zero-energy physical randomness.

When coupled with analog CMOS summing networks, these physical p-bits can solve optimization problems and perform Bayesian inference orders of magnitude faster and with significantly less power than classical CPUs or GPUs.

As we reach the physical limits of traditional semiconductors, shifting from deterministic execution to probabilistic exploitation is no longer just a theoretical curiosity—it is the next frontier of high-performance systems engineering.

#Probabilistic Computing#Neuromorphic Hardware#Algorithms#Python