Simulated Bifurcation (SB) Optimization

• Simulated Bifurcation is a deterministic, quantum-inspired approach that maps combinatorial problems onto a network of nonlinear oscillators using Hamiltonian dynamics.
• It employs a pitchfork bifurcation mechanism and kinetic heating to escape local minima, reducing steps-to-solution in dense Ising and MAX-CUT benchmarks.
• The method leverages hardware accelerations like FPGAs and compute-in-memory architectures, enabling massive parallelism and scalable optimization for large problem instances.

Simulated Bifurcation (SB) is a deterministic, quantum-inspired heuristic for solving large-scale combinatorial optimization problems, specifically those that can be formulated as quadratic unconstrained binary optimization (QUBO) or Ising models. SB operates by simulating the classical Hamiltonian dynamics of a network of coupled nonlinear oscillators whose low-energy phase space encodes optimal or near-optimal solutions to the original combinatorial problem. Its defining features include a pitchfork-type bifurcation mechanism, massive parallelizability, and hardware suitability for high-speed, large-scale integration (Kanao et al., 2022).

1. Mathematical Framework and Core Dynamics

SB targets the minimization of the Ising Hamiltonian

$$E_{\rm Ising} = -\frac{1}{2}\sum_{i=1}^N\sum_{j=1}^N J_{ij} s_i s_j, \qquad s_i\in\{\pm1\},\; J_{ij}=J_{ji},\; J_{ii}=0.$$

Each binary spin $s_i$ is mapped to a pair of continuous dynamical variables: “position” $x_i\in[-1,1]$ and “momentum” $y_i\in\mathbb{R}$ (Kanao et al., 2022). SB dynamics are governed by Hamiltonian-like ODEs, with two canonical variants: $$\begin{cases} \dot{x}_i = a_0 y_i, \ \dot{y}_i = -[a_0 - a(t)] x_i + c_0 f_i, \end{cases}$$ where

$$f_i = \begin{cases} \sum_{j=1}^N J_{ij} x_j & \text{(ballistic SB, bSB)} \ \sum_{j=1}^N J_{ij} \operatorname{sgn}(x_j) & \text{(discrete SB, dSB)} \end{cases}$$

The bifurcation parameter $a(t)$ is ramped from $0$ to $a_0$ over the simulation (Kanao et al., 2022). State variables are confined to “inelastic walls” $|x_i|\le1$, and the discrete solution is recovered by $s_i = \operatorname{sgn}(x_i)$ at termination.

Generalization to higher-order cost functions (e.g., cubic terms) is achieved by embedding the corresponding polynomial gradient directly in $f_i$ (Kanao et al., 2022). For constraint-laden QUBO mappings (e.g., community detection, traffic signal control), the Ising mapping can become high-dimensional, but the SB update equations remain structurally consistent (Li et al., 2024, Zhao et al., 18 Feb 2025).

2. Extensions: Thermalization and the Heated SB Algorithm

A deterministic emulation of thermal noise is achieved by augmenting the equations with a Nosé–Hoover thermostat: $$\begin{cases} \dot{x}_i = a_0 y_i, \ \dot{y}_i = -[a_0 - a(t)] x_i + c_0 f_i - \xi y_i, \ \dot{\xi} = \frac{1}{M}\left(\sum_{i=1}^N y_i^2 - N T\right), \end{cases}$$ where $T$ is the target temperature (Kanao et al., 2022).

Empirically, only the heating phase ($T_{\mathrm{inst}} < T$) enhances exploration; for large $M$, $\xi$ is nearly constant and can be replaced by $-\gamma$ ($\gamma>0$), yielding the Heated Simulated Bifurcation (HSB) equations: $$\dot{y}_i = -[a_0 - a(t)] x_i + c_0 f_i + \gamma y_i.$$ The $+\gamma y_i$ term injects kinetic energy (deterministically) into each oscillator, enabling escape from local minima in the Ising energy landscape. HSB achieves higher maximum and average cut values than unheated SB, with up to 33% reduction in steps-to-solution for $N=2000$ fully connected Ising problems (Kanao et al., 2022).

3. Algorithmic Schemes, Complexity, and Benchmarks

Algorithmic steps for SB (e.g., symplectic Euler integration, “perfectly inelastic” boundary enforcement) are explicitly amenable to high parallelism (Li et al., 2024, Orlando et al., 14 Oct 2025, Kashimata et al., 2023). Key computational features:

• Per-step complexity: $O(N^2)$ for dense Ising spin networks due to matrix-vector products.
• For problems with reducible coupling structure $J=WW^\top$ and $W\in\mathbb{R}^{N \times M}$, the “high-dimensional expansion” lowers the step cost to $O(N^{1.35})$ with $M=O(N^{0.35})$ (Zhao et al., 18 Feb 2025).
• Hardware implementations with fixed-point (FPGA) or compute-in-memory (CIM) architectures further optimize the update loop for resource and energy constraints (Orlando et al., 14 Oct 2025, Dee et al., 14 Apr 2025, Qian et al., 19 Dec 2025).
• Solution performance for dense MAX-CUT instances with $N=2000$ (Sherrington–Kirkpatrick model): HSB achieves a probability $P=0.47$ of finding the optimum at $N_s=10^4$ steps, with $32.7\%$ fewer steps-to-solution compared to dSB; mean cut value $C=33250$ vs. $33190$ for dSB (Kanao et al., 2022).

Table: Empirical SB/HSB Performance, $N=2000$ MAX-CUT (Kanao et al., 2022)

Algorithm	Avg $C$	Max $C$	$P$ at $N_s=10^4$	Step-to-solution $S$
bSB	33050	33270	0.00	—
dSB	33190	33337	0.30	$1.23 \times 10^4$
HbSB	33250	33337	0.47	$8.27 \times 10^3$
HdSB	33230	33337	0.38	$9.88 \times 10^3$

4. Dynamical Properties and Generalized Controls

The success of SB hinges on controlled dynamical evolution and the pitchfork bifurcation induced by the time-varying parameter $a(t)$ (or, in generalized formulations, $p(t)$). The original “ballistic” SB uses a global $a(t)$, but generalized SB (GSB) introduces per-oscillator nonlinear control parameters $p_i(t)$: $$p_i(t_{m+1}) = p_i(t_m) - [1 - A x_i^2(t_m)] \frac{p_i(t_m)}{M - m},$$ where $A$ is the nonlinear-control strength. $A=0$ recovers standard SB; $A>0$ adaptively slows $p_i$’s schedule for oscillators “stuck” at boundaries, thus assisting escapes from shallow local minima (Goto et al., 25 Aug 2025).

GSB can reach $P_S \approx 100\%$ on $N=2000$ random dense graphs, achieving time-to-solution of 9.6 ms (FPGA), two orders of magnitude faster than SB (1.3 s) (Goto et al., 25 Aug 2025). Optimal $A$ operates at the “edge of chaos,” where trajectory sensitivity is maximized but not fully randomized.

5. Hardware Realizations and Parallel Implementations

SB’s suitability for hardware stems from its reliance on dense matrix–vector products, thresholding, and simple nonlinearities. Implementation developments include:

• Discrete SB (dSB) and its heating variant (HdSB), which use sign-based couplings, are highly compatible with fixed-point and integer hardware (e.g., FPGAs, SRAM/FeFET CIM) (Orlando et al., 14 Oct 2025, Dee et al., 14 Apr 2025, Qian et al., 19 Dec 2025).
• Mixed-signal SRAM-CIM solvers achieve sub-microsecond (0.6 μs, 60 spins) time-to-solution at $\sim$10 mW power, surpassing previous digital and analog Ising solvers (Dee et al., 14 Apr 2025).
• Multi-chip SB machines with dual-ring streaming architectures exhibit near-ideal strong scaling ($97.9\%$ pipeline efficiency) for up to 100,000 spins, matching state-of-the-art optical Ising performance (Kashimata et al., 2023).
• Recent BEOL FeFET-CIM SB solvers integrate attention-inspired initialization and lightweight update rules (“Light SB”) for ultra-fast and energy-efficient optimization at scale (Qian et al., 19 Dec 2025).

6. Application Domains

SB has been effectively applied across a broad range of combinatorial optimization problems:

• MAX-CUT / Ising benchmarks: Dense and sparse graph variants, including Sherrington–Kirkpatrick and G-set (Kanao et al., 2022, Goto et al., 25 Aug 2025).
• Community detection: QUBO-reformulated modularity maximization; empirical results match or exceed Fujitsu’s Digital Annealer and outperform D-Wave/IBM quantum solvers (Li et al., 2024).
• MIMO detection: SB achieves significantly lower BER than MMSE and coherent Ising machine (CIM) methods and can be further enhanced by deep-unfolded (differentiable) variants and Levenberg–Marquardt preconditioning (Zhang et al., 2022, Takabe, 2023).
• Portfolio optimization: Efficient solutions to large-scale integer-constrained Markowitz problems and trading trajectory optimization, scaling to $N\sim 1000$ assets in subsecond wall-time (Bouquet et al., 2021, Steinhauer et al., 2020).
• Traffic signal control: Real-world networks up to $N=604$ intersections; high-dimensional expansion reduces complexity to $O(N^{1.35})$, outperforming Metropolis SA in both runtime and final objective (Zhao et al., 18 Feb 2025).
• Higher-order Ising terms: Third-order SB (3bSB) handles sparse polyadic constraints directly, outperforming reduction-to-second-order+gadgets and simulated annealing (Kanao et al., 2022).

7. Limitations, Trade-offs, and Prospects

SB’s deterministic, fully parallel dynamics lead to several salient features:

• Deterministic runs, unlike simulated annealing’s (SA) stochasticity.
• Absence of random numbers in HSB: momentum “heating” is sufficient for explorative search (Kanao et al., 2022).
• Wall boundary enforcement prevents variable escape but can induce trapping in absence of “jitter.”
• Parameter sensitivities: integrator step size ($\Delta t$), bifurcation schedule ($a(t)$), and heating rate ($\gamma$), require tuning per problem class.
• The generalized and hardware-enabled SB family continues to evolve, with tunable chaos (GSB) and initialization heuristics (attention-inspired, CiM-friendly) to further shorten convergence time and improve robustness (Goto et al., 25 Aug 2025, Qian et al., 19 Dec 2025).

A plausible implication is that further extensions (adaptive schedules, hybrid stochastic-deterministic control, tunable network structure) and continued hardware–algorithm co-design will broaden SB’s applicability to harder, more structured classes of NP-hard combinatorial problems at industrial scales (Goto et al., 25 Aug 2025, Kashimata et al., 2023, Qian et al., 19 Dec 2025).