Simulating Quantum Decision Models in Pure Python (No QPU Required)

I spent most of last week arguing with a vendor who insisted we needed cloud QPU access to run our new multi-agent decision matrix. Actually, let me back up — I told them to shove it.

Look, the hype around quantum hardware right now is exhausting. But you know, everyone assumes that if you want to use quantum mechanics to model complex, overlapping probabilities—like evaluating a robot’s actions for safety, fairness, and efficiency simultaneously—you need a physical quantum computer. You don’t. You just need the math.

Quantum mechanics, at its software layer, is just linear algebra. It’s complex numbers and matrices. And you can run that entirely classically in pure Python. I recently ported a massive “quantum-inspired” ethical evaluator down to raw NumPy. It evaluates high-dimensional decision spaces without ever touching a noisy, expensive quantum processor.

And here is exactly how I built it, and why you should probably cancel your quantum cloud subscription if you’re just doing state-space modeling.

The Math Behind the Magic

When we talk about a “quantum” approach to decision-making, we’re usually talking about superposition and interference. Instead of a classical decision tree where an AI evaluates path A, then path B, then path C, we represent all possible decisions as a single state vector.

We then apply unitary matrices (operators) that act as our constraints. Safe actions get constructive interference (their probabilities amplify). Risky actions get destructive interference (they cancel out). Then we “measure” the system to get the optimal choice.

You can simulate this perfectly on a MacBook. The only catch is memory. A quantum system with $N$ variables requires an array of size $2^N$. For 50 variables, you’d need a supercomputer. But for a highly targeted 12-factor alignment core? That’s an array of 4,096 complex floats. Your phone could calculate that in milliseconds.

Building the Classical-Quantum Core

I tested this on Python 3.12.2 using numpy 1.26.4. The goal was to build a pure Python class that initializes a decision space, throws it into superposition, applies risk penalties via phase shifts, and collapses to the safest action.

import numpy as np

class ClassicalQuantumCore:
    def __init__(self, num_factors):
        """
        Initialize a 2^N dimensional Hilbert space.
        num_factors = the number of binary decisions/constraints.
        """
        self.n = num_factors
        self.dim = 2**self.n
        
        # Start in the |0...0> state
        self.state = np.zeros(self.dim, dtype=np.complex128)
        self.state[0] = 1.0 + 0.0j

    def apply_superposition(self):
        """
        Applies a Hadamard transform to put all possible 
        decision states into equal superposition.
        """
        # The basic 1-qubit Hadamard matrix
        H = np.array([[1, 1], [1, -1]]) / np.sqrt(2)
        
        # Tensor product it N times to cover the whole space
        transform = H
        for _ in range(self.n - 1):
            transform = np.kron(transform, H)
            
        # Matrix multiplication to update the state
        self.state = transform @ self.state

    def apply_risk_penalty(self, risk_profile):
        """
        risk_profile: A 1D array of length 2^N containing risk scores (0.0 to 1.0)
        We convert these to phase shifts. High risk = destructive interference.
        """
        # Convert risk to an angle between 0 and Pi
        phases = np.exp(1j * risk_profile * np.pi)
        
        # Apply the phase shift (diagonal operator)
        self.state = self.state * phases

    def measure(self):
        """
        Collapse the state vector into classical probabilities.
        """
        probabilities = np.abs(self.state)**2
        
        # Return the index of the highest probability state
        # (This represents the optimal, safest decision path)
        return np.argmax(probabilities), probabilities

This is obviously a simplified version of the math, but the skeleton is exactly what runs in production. You initialize the space. You spread the probability across all possible actions. You apply your constraints as phase shifts. You measure the result.

A Massive Gotcha: The np.kron Trap

And I need to point out something that cost me an entire afternoon. In the code above, I used np.kron (Kronecker product) in a loop to build the superposition transform. Do not do this in production if $N > 10$.

If you have 12 factors, np.kron builds a $4096 \times 4096$ dense matrix before doing the multiplication. I crashed my M3 Max running out of memory when I tried pushing it to 18 factors. The naive matrix multiplication scales horribly.

The fix? You don’t actually need to build the giant matrix. You can apply the 2×2 Hadamard matrix to each “qubit” sequentially by reshaping the state vector. Here’s the optimized method I ended up writing:

    def apply_superposition_fast(self):
        """
        Memory-efficient Hadamard application without building 
        the massive 2^N x 2^N operator matrix.
        """
        H = np.array([[1, 1], [1, -1]]) / np.sqrt(2)
        
        # Reshape state to a multi-dimensional tensor
        tensor_state = self.state.reshape([2] * self.n)
        
        # Apply H to each dimension using einsum
        for i in range(self.n):
            # Move the target dimension to the front, multiply, put it back
            tensor_state = np.tensordot(H, tensor_state, axes=(1, i))
            tensor_state = np.moveaxis(tensor_state