Coupled-Field Dynamical Relaxation for QUBO and Ising Optimizations

1. Introduction

Quadratic Unconstrained Binary Optimization (QUBO) problems arise in logistics, medical imaging, portfolio optimization, robotics, scheduling, and machine learning [1,2,3,4,5]. The standard QUBO objective seeks to minimize $f\left(x\right) = x^T Q x$ over binary variables x ∈ {0,1}^N, which can be equivalently expressed as an Ising Hamiltonian $H = {\sum }_{ij}{J}_{ij} s_i s_j+ {\sum }_i{h}_i s_i$ with spins s_i ∈ {-1,+1}. Since general QUBO/Ising optimization is NP-hard, practical systems rely on specialized heuristics or hardware accelerators.

Quantum annealers (D-Wave) and coherent optical Ising machines implement energy minimization through quantum dynamics or nonlinear photonics [6,7,8,9,10], but require cryogenic conditions or ultrastable lasers, limiting accessibility and scalability. This work introduces a fully classical, deterministic architecture operating at room temperature with no quantum coherence requirement, yet retaining the essential property needed for optimization: a multi-element system that naturally evolves toward energy minima.

1.1. Novel Contributions and Scope

This work contributes a theoretical framework with four distinguishing features not present in existing oscillator-based Ising machines:

Unified Lyapunov Energy Functional: CF-bit dynamics are derived from a single explicit energy functional $E_{total} = \sum {\sum }_{i}V\left({\phi }_i\right) + {\sum }_{ij}{J}_{ij} cos\left({\phi }_i- {\phi }_j\right)$, enabling analytical convergence guarantees. Existing optical and electronic oscillator machines do not provide a closed-form global Lyapunov function.

Binary Encoding in the Potential Itself: The binary state arises from minima of the intrinsic potential $V\left(\phi \right) = -cos\left(2\phi \right)$, not from external thresholding or gain saturation. This enables a clean algebraic mapping from QUBO coefficients to coupling strengths.

Deterministic Lyapunov Dynamics: The CF system guarantees monotonic energy descent dE/dt ≤ 0 under damping, unlike coherent Ising machines where gain saturation and noise injection are intrinsic. Stochastic noise is not required for stability or final convergence but may optionally enhance exploration.

Direct Algebraic QUBO Mapping: The transformation from QUBO matrix Q to CF couplings {J_ij, h_i} is derived analytically from the energy functional, not heuristically fit as in some optical implementations.

Scope and Limitations: This work is presented as a theoretical and conceptual contribution emphasizing dynamical structure and physical interpretability. Numerical demonstrations are limited to instances with N ≤ 20 bits. The framework does not overcome NP-hardness and may converge to local minima for frustrated systems. Larger-scale behavior, physical coupling topology costs, and competitive benchmarking against state-of-the-art solvers remain important directions for future research but lie beyond the scope of the present study.

2. Coupled-Fields Theory and the CF-Bit

2.1. Derivation from Energy Functional

A CF-bit is a classical nonlinear oscillator with internal phase state φ_i evolving in a bistable periodic potential. We begin from a Lyapunov-like energy functional for a single CF-bit:

$$E_{single} = V\left(\phi \right) + (m/2){(d\phi /dt)}^2$$

(1)

where V(φ) = -cos(2φ) defines the bistable potential and m represents effective inertia inherited from underlying field dynamics. The equation of motion follows from the Euler-Lagrange formulation with linear damping:

$$m d^2\phi /{dt}^2 + \alpha d\phi /dt = -dV/d\phi$$

(2)

Substituting V(φ) = -cos(2φ) yields:

$$m d^2\phi /{dt}^2 + \alpha d\phi /dt = -2 sin\left(2\phi \right)$$

(3)

This second-order differential equation describes a damped nonlinear oscillator. For the binary case (n=2 in the general formulation $V\left(\phi \right) = -cos\left(n\phi \right))$, the potential exhibits exactly two stable minima at φ = 0 and φ = π, corresponding to logical states 0 and 1 (see Figure 1).

2.2. Network-Level Lyapunov Stability

For a network of N coupled CF-bits, the total energy functional is:

$$E_{total} = {\Sigma }_i V\left({\phi }_i\right) + \left({\displaystyle \frac{m}{2}}\right){\Sigma }_i {(d{\phi }_i/dt)}^2 + {\sum }_{i<j} J_{ij}cos\left({\phi }_i- {\phi }_j\right) + {\Sigma }_i{h}_i{\phi }_i$$

(4)

This serves as a global Lyapunov function. Computing its time derivative:

$$dE/dt = {\Sigma }_i \left[\right(dV/d{\phi }_i + \partial E_{coupling}/\partial {\phi }_i\left)\right(d{\phi }_i/dt) + m(d^2{\phi }_i/{dt}^2\left)\right(d{\phi }_i/dt\left)\right]$$

(5)

Using the equation of motion

$$m d^2{\phi }_i/{dt}^2=-(dV/d{\phi }_i)-(\partial E_{coupling}/\partial {\phi }_i)-\alpha (d{\phi }_i/dt)$$

(6)

we obtain:

dE/dt = -α Σ_i (dφ_i/dt)² ≤ 0

(7)

This confirms monotonic energy descent for α > 0, guaranteeing convergence to fixed points. This Lyapunov property distinguishes the CF architecture from coherent Ising machines (CIMs) and oscillator-based Ising machines (OIMs), where dynamics include non-gradient terms (gain saturation, noise) and no global Lyapunov functional has been established [11,12,13].

3. Interaction Between CF-Bits and Ising Mapping

CF-bits interact through coupling terms depending on phase differences:

$$E_{coupling} = {\sum }_{i<j} J_{ij}cos\left({\phi }_i- {\phi }_j\right)$$

(8)

Near the stable phases φ_i ∈ {0, π}, expanding to second order around these minima and defining binary variables $b_i$ via ${\phi }_i \approx b_i \pi$(with $b_i \in \{0,1\}$), the coupling energy reduces to:

$$E_{coupling} \approx {\Sigma }_{i<j} J_{ij} cos(b_i \pi - b_j \pi ) = {\Sigma }_{i<j}{J_{ij}\left(-1\right)}^{(b_i + b_j)}$$

(9)

Using the Ising transformation s_i = 2b_i - 1 (mapping {0,1} to {-1,+1}), this becomes:

$$E_{Ising} = {\Sigma }_{i<j}{s}_i s_j + {\Sigma }_i{h}_i s_i \left(9\right)$$

(10)

demonstrating that the CF network energy naturally reduces to an Ising Hamiltonian. Local fields h_i can be encoded through bias terms added to individual CF-bit potentials. This provides the formal justification for mapping QUBO coefficients directly to CF coupling strengths.

4. Mapping QUBO to CF Networks

A QUBO problem seeks to minimize $f\left(x\right) = x^T Q x$ with $x \in {\left\{0,1\right\}}^N$. Using the transformation $s_i = 2x_i - 1$, the QUBO objective becomes equivalent to the Ising Hamiltonian:

$$H_{Ising} = \left({\displaystyle \frac{1}{4}}\right){\Sigma }_i {\Sigma }_j Q_{ij} s_i s_j + \left({\displaystyle \frac{1}{2}}\right){\Sigma }_i \left({\Sigma }_j Q_{ij}\right)s_i + constant \left(10\right)$$

(11)

Identifying the coupling matrix and local fields:

$$\begin{gathered} J_{ij} ={\displaystyle \frac{Q_{ij}}{4}} \left(i \ne j\right) \\ {h}_i = \left({\displaystyle \frac{1}{2}}\right){\Sigma }_j Q_{ij} \end{gathered}$$

(12)

This mapping is lossless: every QUBO instance produces a unique CF coupling matrix whose minima correspond exactly to QUBO minima. Since each CF-bit has exactly two stable attractors (φ = 0, π), the transformation guarantees one-to-one correspondence between QUBO solutions and CF equilibrium configurations.

4.1. Explicit 3×3 Example

To illustrate the mapping concretely, consider the 3×3 QUBO matrix Q:

$$Q = \left[\begin{array}{ccc}2& -1& 0\\ -1& 3& -2\\ 0& -2& 1\end{array}\right]$$

(13)

The QUBO cost function is $f\left(x\right) = x^T Q x$. Applying the binary-to-Ising transformation s_i = 2x_i - 1 and expanding yields an Ising energy with couplings and fields determined analytically from Q entries. For this specific matrix:

$$\begin{gathered} J_{12} = -1/4, J_{13} = 0, J_{23} = -1/2 \\ {h}_1 = 1/2, h_2 = 0, h_3 = -1/2 \end{gathered}$$

These Ising parameters map directly to the CF network energy functional, defining the physical couplings and bias coefficients. The minima of E_CF are in one-to-one correspondence with the minima of the original QUBO function f(x).

5. Dynamics and Convergence

The time evolution of each CF-bit follows:

$$m d^2{\phi }_i/{dt}^2 + \alpha d{\phi }_i/dt = -2sin\left(2{\phi }_i\right) + {\Sigma }_j J_{ij} sin\left({\phi }_i- {\phi }_j\right) - h_i \left(14\right)$$

(14)

For numerical simulations, the network was integrated using a fourth-order Runge-Kutta scheme with fixed time step dt = 0.01. The damping parameter α = 1.1 was chosen to ensure monotonic energy decrease while allowing sufficient transient oscillation to traverse shallow minima. Convergence was defined as $|d{\phi }_i/dt| < \epsilon$ for all i, with $\epsilon = {10}^{-3}$.

Key dynamical features:

Parallel evolution: All CF-bits update simultaneously, providing massive analog parallelism.

Energy descent: Damping guarantees dE/dt ≤ 0, ensuring relaxation toward energy minima.

Deterministic operation: Noise is not intrinsic to stability or late-time phase locking. However, optional transient stochastic perturbations can enhance exploration of rugged landscapes.

Multi-minima exploration: Large initial oscillations enable the system to bypass shallow minima during early evolution.

6. Simulation Results

We present numerical demonstrations on three problem classes. All simulations used the parameters listed in Table 1 below.

We present numerical demonstrations on three problem classes. All simulations used the following parameters unless otherwise noted:

Integration method: 4th-order Runge-Kutta

Time step: dt = 0.01

Damping coefficient: α = 1.1

Initial phases: Uniform random on [0, 2π]

Initial velocities: Gaussian with σ = 0.05

Convergence threshold: $\epsilon = {10}^{-3}$ QUBO coefficient range: Normalized to [-1, 1]

6.1. 3-Bit Illustrative Example

A 3-bit QUBO with Q given in Equation (12) was encoded into a 3-CF-bit system. The system evolved deterministically from random initial conditions to a final configuration φ_final ≈ (π, π, 0), corresponding to binary solution x = (1, 1, 0). This matches one of the global minima found by brute-force enumeration, validating the mapping.

Figure 4 shows the energy convergence (top panel) and phase trajectories (bottom panel). The CF energy decreases monotonically from its initial value to a stable minimum within approximately 150 iterations, demonstrating the Lyapunov descent property.

Figure 2. Conceptual schematic of a single CF-bit oscillator. The CF-bit consists of an LC resonator with intrinsic binary phase states, a damping resistor R enforcing relaxation toward stable equilibria, and optional input for initialization or coupling. This illustrates one possible classical implementation pathway.

Figure 2. Conceptual schematic of a single CF-bit oscillator. The CF-bit consists of an LC resonator with intrinsic binary phase states, a damping resistor R enforcing relaxation toward stable equilibria, and optional input for initialization or coupling. This illustrates one possible classical implementation pathway.

Figure 3. Mapping a QUBO matrix onto a network of interacting CF-bits. Each binary variable x_i is encoded as a CF-bit with internal phase φ_i. Quadratic coefficients Q_ij translate into pairwise coupling strengths J_ij, while linear terms become local bias fields h_i.

Figure 3. Mapping a QUBO matrix onto a network of interacting CF-bits. Each binary variable x_i is encoded as a CF-bit with internal phase φ_i. Quadratic coefficients Q_ij translate into pairwise coupling strengths J_ij, while linear terms become local bias fields h_i.

Figure 4. CF convergence on a 3-bit QUBO instance. Top panel: CF energy decreases monotonically as damped dynamics relax toward a minimum. Bottom panel: Phase trajectories φ_i(t) starting from random initial conditions and settling near 0 or π, corresponding to discrete bit states encoding a QUBO minimum.

Figure 4. CF convergence on a 3-bit QUBO instance. Top panel: CF energy decreases monotonically as damped dynamics relax toward a minimum. Bottom panel: Phase trajectories φ_i(t) starting from random initial conditions and settling near 0 or π, corresponding to discrete bit states encoding a QUBO minimum.

6.2. Random 20-Bit QUBO Instances

To probe larger systems, twenty dense 20-bit QUBO instances were generated with random symmetric matrices Q. For each instance, the corresponding Ising parameters were computed and mapped to a 20-CF-bit network. The system was integrated until convergence or until a maximum of 500 iterations.

Results are summarized in Table 2. The CF network reached the global optimum (verified by brute-force enumeration) in 15 out of 20 cases (75% success rate). The remaining 5 instances converged to local minima one energy level above the global minimum. Mean convergence time was 171 iterations (σ = 22 iterations).

This performance is consistent with other analog Ising-type solvers and demonstrates that CF dynamics provide an effective heuristic for mid-size QUBO instances. The deterministic relaxation properties ensure stable convergence, while the continuous phase space allows exploration of the energy landscape during transient evolution.

Summary: Global optimum reached in 15/20 instances (75%). Near-optimal solutions (within one energy level) in the remaining 5/20 (25%). Mean convergence: 171 iterations. This behavior is consistent with analog Ising-type dynamics exhibiting monotonic relaxation with occasional local minima.

6.3. Structured MaxCut Instances

Structured QUBO instances from MaxCut problems on graphs with 10-20 vertices were tested. For these problems, the CF network exhibited rapid initial energy descent followed by stable phase-locking. Final bit assignments corresponded to cuts matching optimal MaxCut values (for small graphs verified by exhaustive search) or high-quality solutions comparable to established heuristics for larger graphs.

Increasing initial oscillation amplitude or adjusting the damping parameter improved escape from shallow local minima, analogous to temperature schedule tuning in simulated annealing. These structured problems demonstrate that CF dynamics successfully reproduce behavior expected from physical Ising hardware on both random and structured optimization instances.

6.4. Computational Cost and Scaling

Each integration step requires O(N²) operations for dense coupling matrices (evaluating all pairwise interactions) or O(N·k) for sparse instances with k neighbors per CF-bit. For the 20-bit problems tested, typical convergence required 150-200 iterations, giving total computational cost of approximately 30,000-40,000 pairwise evaluations.

These costs reflect numerical integration of continuous dynamics rather than optimized algorithmic implementations. For comparison, Table 3 shows solution quality and approximate runtime for different approaches on the same 20-bit test set. Note that these comparisons are provided solely as qualitative reference for small problem instances and should not be interpreted as competitive benchmarks against mature solvers.

The CF approach achieves solution quality comparable to simulated annealing with faster convergence (measured in integration steps), though this comparison reflects research code implementations rather than optimized production solvers. Systematic scaling analysis beyond N=20 and comparison with state-of-the-art optimization libraries remain important directions for future work.

7. Comparison with Alternative Optimization Approaches

7.1. Relation to Existing Oscillator-Based Ising Machines

Several prior approaches use nonlinear phase dynamics for optimization, including LC oscillator networks [11], MEMS-based Ising solvers, and coherent optical machines [8,9,10]. The CF architecture differs in three key respects:

CF dynamics are derived from an explicit global Lyapunov functional E_total, enabling analytical convergence guarantees and transparent energy-descent interpretation.

The CF potential directly embeds binary states through intrinsic bistability V(φ) = -cos(2φ), yielding a clean algebraic mapping to the Ising Hamiltonian without external thresholding.

The architecture operates deterministically without requiring gain saturation, stochastic noise injection, or optical coherence for stability or final convergence.

These distinctions clarify the theoretical contribution: the CF model provides a classical, room-temperature framework with mathematically transparent dynamics and direct QUBO-to-coupling mapping.

7.2. Non-Existence of Lyapunov Functional in CIM/OIM Architectures

A central distinction between the CF architecture and coherent Ising machines (CIMs) or oscillator-based Ising machines (OIMs) is the existence of a global Lyapunov functional governing the dynamics.

CIM dynamics based on degenerate optical parametric oscillators (DOPOs) evolve according to equations that include gain-saturation and pump-depletion terms [8]. These terms cannot be expressed as gradients of any global potential function E. As noted in [8], CIM dynamics are fundamentally non-gradient, relying on stochastic amplitude bifurcation to reach low-energy states. Similarly, OIM systems based on coupled LC oscillators include phase-normalization dynamics and amplitude-regulating terms that violate gradient structure [11,12].

The consequence is that CIM and OIM trajectories may temporarily increase in effective Ising energy, and convergence relies on noise, bifurcation, or gain-control mechanisms rather than guaranteed monotonic relaxation. In contrast, the CF architecture possesses an explicit Lyapunov functional whose gradient fully determines the dynamics, ensuring dE/dt ≤ 0 at all times under damping.

This structural property—to our knowledge not demonstrated for optical or electronic oscillator Ising machines in their standard formulations—provides the CF framework with mathematically transparent convergence behavior and enables design principles based on energy-landscape engineering.

7.3. Relation to Classical Neural Energy Models

Early neural energy models such as Hopfield networks [14] and Cohen-Grossberg systems [15] demonstrated that gradient-based dynamics can perform combinatorial optimization by descending a Lyapunov energy landscape. The CF model shares this gradient-descent principle but differs in implementation:

Hopfield networks use threshold units with discrete states, while CF-bits are continuous phase oscillators that naturally explore the energy landscape during transient evolution.

The CF architecture derives binary computation from intrinsic potential bistability V(φ) = -cos(2φ), whereas Hopfield networks impose discrete thresholds externally.

CF dynamics include second-order (inertial) terms, enabling richer transient behavior and natural frequency-based coupling mechanisms.

The CF framework can thus be viewed as a classical oscillator-based realization of energy-descent optimization, bridging concepts from neural computation and analog physical systems.

8. Convergence Properties and Role of Noise

8.1. Deterministic Lyapunov Convergence

Under purely deterministic damped dynamics (no stochastic perturbations), the CF system exhibits monotonic energy descent $dE/dt = -\alpha {\Sigma }_i {(d{\phi }_i/dt)}^2 \le 0$. This guarantees convergence to a fixed point corresponding to a local (or global) minimum of the energy landscape. However, as with all gradient-descent methods on non-convex landscapes, deterministic dynamics may become trapped in the nearest local minimum.

8.2. Noise-Assisted Exploration

To enhance exploration of rugged energy landscapes, optional transient stochastic perturbations can be introduced:

$$m d^2{\phi }_i/{dt}^2 + \alpha d{\phi }_i/dt = -dE/d{\phi }_i + \sigma \left(t\right)·{\xi }_i\left(t\right)$$

(15)

where ${\xi }_i\left(t\right)$ is white Gaussian noise and σ(t) is a time-dependent noise amplitude. The noise amplitude is annealed according to $\sigma \left(t\right) = {\sigma }_0 exp(-t/\tau )$, starting at ${\sigma }_0 \approx 0.1-0.5$ and decaying with time constant τ chosen to match the typical relaxation time of the system.

During the exploratory regime (large σ), the noise temporarily breaks Lyapunov monotonicity, enabling barrier crossing and escape from shallow local minima. As σ → 0 in the late-time regime, the system transitions back to deterministic convergence, ensuring stable final phase-locking. This approach combines the stability guarantees of deterministic Lyapunov dynamics with enhanced global search capability.

Figure 6 illustrates this behavior: the deterministic case (smooth curve) shows strict energy monotonicity, while the noise-assisted case (fluctuating curve) exhibits temporary energy increases during exploration followed by convergence to a deeper minimum as noise is reduced.

Figure 5. CF convergence on a 20-bit QUBO instance showing energy vs. time. The smooth monotonic descent reflects the Lyapunov property dE/dt ≤ 0, with eventual plateau indicating convergence to a stable low-energy state.

Figure 5. CF convergence on a 20-bit QUBO instance showing energy vs. time. The smooth monotonic descent reflects the Lyapunov property dE/dt ≤ 0, with eventual plateau indicating convergence to a stable low-energy state.

Figure 6. Effect of noise on convergence dynamics. Deterministic case (smooth curve) shows strict monotonic descent. Noise-assisted case (fluctuating curve) exhibits temporary energy increases enabling barrier crossing, followed by convergence to a deeper minimum as noise amplitude decays.

Figure 6. Effect of noise on convergence dynamics. Deterministic case (smooth curve) shows strict monotonic descent. Noise-assisted case (fluctuating curve) exhibits temporary energy increases enabling barrier crossing, followed by convergence to a deeper minimum as noise amplitude decays.

8.3. Comparison with Coherent Ising Machines

In coherent Ising machines (CIMs), stochastic noise from quantum vacuum fluctuations and amplifier noise is intrinsic to the system dynamics and required for optimization performance [8]. The CIM phase dynamics can be written schematically as dφ_i/dt = F_i({φ_j}) + ξ_i(t), where ξ_i represents unavoidable stochastic processes. No global Lyapunov functional exists, and convergence is fundamentally probabilistic.

In contrast, the CF model separates deterministic convergence (guaranteed by the Lyapunov functional) from optional noise-assisted exploration (controllable and transient). This design choice enables stable late-time phase locking while retaining the ability to enhance global search when needed. The noise amplitude and annealing schedule are external control parameters rather than intrinsic physical constraints.

9. Applications

The CF architecture is applicable to combinatorial optimization problems that can be formulated as QUBO or Ising instances, including:

Logistics and routing (vehicle routing, traveling salesman)

Portfolio optimization and financial risk minimization

Scheduling and resource allocation

Graph problems (MaxCut, graph coloring, maximum independent set)

Machine learning (training binary neural networks, feature selection)

VLSI design and circuit layout optimization

The specific advantage of the CF approach lies in its deterministic convergence guarantees, room-temperature operation, and transparent energy-landscape structure, which may facilitate analysis and design of problem-specific coupling architectures.

10. Discussion and Future Directions

This work introduces a coupled-field dynamical framework in which discrete optimization emerges from continuous relaxation of phase fields governed by an explicit global energy functional. The key theoretical contributions are (i) derivation of dynamics from a single Lyapunov function enabling analytical convergence analysis, (ii) natural algebraic mapping from QUBO to CF couplings, and (iii) deterministic operation without intrinsic noise requirements.

Important directions for future research include:

Systematic scaling analysis up to N=100-1000 bits to characterize performance degradation and identify problem classes where CF dynamics remain effective.

Computational complexity characterization: theoretical analysis of iteration count vs. problem size and coupling density.

Physical implementation feasibility studies examining electronic, mechanical, or photonic oscillator realizations and their engineering constraints.

Hybrid classical-CF algorithms combining CF relaxation with discrete search heuristics.

Analysis of coupling topology costs for dense vs. sparse QUBO matrices and development of approximation schemes for hardware-constrained connectivity.

Application-specific architectures optimized for particular problem classes (e.g., MaxCut, SAT, scheduling).

10.1. Limitations

The CF architecture does not overcome NP-hardness and may converge to local minima for frustrated systems with conflicting constraints. Numerical demonstrations are limited to N ≤ 20 bits; behavior at larger scales requires further investigation. Physical implementation of fully dense coupling matrices scales as O(N²) in component count, potentially limiting scalability without sparse approximations or hierarchical architectures.

These limitations are characteristic of analog Ising-type solvers and define important trade-offs between problem size, coupling density, solution quality, and hardware complexity.

11. Conclusions

We have introduced a coupled-field dynamical framework in which discrete optimization problems are represented as energy landscapes governing the relaxation of interacting phase fields. By constructing an explicit global energy functional and deriving overdamped gradient dynamics, the framework admits a Lyapunov description guaranteeing monotonic energy descent and stable late-time phase locking.

This structure provides a physically transparent alternative to oscillator-based Ising approaches whose dynamics include non-gradient terms and for which no global Lyapunov functional has been established. The analysis clarifies that stochastic perturbations are not intrinsic to stability or final convergence but can be introduced transiently to enhance exploration of complex landscapes.

Numerical results on instances up to 20 bits demonstrate the dynamical mechanism and validate the QUBO-to-CF mapping. While the present study is theoretical in nature and limited in scale, it highlights how discrete computational behavior can arise from continuous nonlinear dynamics governed by a single energy functional. Extensions to larger problem instances, systematic performance studies, and physical realizations constitute important directions for future work.