The Simulated Bifurcation (SB) algorithm is a fast and highly parallelizable state-of-the-art algorithm for combinatorial optimization inspired by quantum physics and spins dynamics. It relies on Hamiltonian quantum mechanics to find local minima of Ising problems. The last accuracy tests showed a median optimality gap of less than 1% on high-dimensional instances.
This open-source package utilizes PyTorch to leverage GPU computations, harnessing the high potential for parallelization offered by the SB algorithm.
It also provides an API to define Ising models or other NP-hard and NP-complete problems (QUBO, Karp problems, ...) that can be solved using the SB algorithm.
Compute Plateform | CPU | GPU |
---|---|---|
Instructions |
pip install simulated-bifurcation |
Install PyTorch with GPU support pip install simulated-bifurcation |
An Ising problem, given a null-diagonal square symmetrical matrix
This problem is known to be NP-hard but is very useful since it can be used in many sectors such as finance, transportation or chemistry or derived as other well-know optimization problems (QUBO, MAXCUT, Knapsack problem, etc.).
The Simulated Bifurcation algorithm was originally introduced to solve Ising problems by simulating the adiabatic evolution of spins in a quantum Hamiltonian system, but can also be generalized to a wider range of optimization problems.
The SB algorithm can be written as the minimization or maximization of multivariable polynomials of degree two, i.e. written as
for which the
This can also be seen as the sum of a quadratic form, a linear form and a constant term and such a formulation is the basis of many optimization problems.
The minimize
and maximize
functions allow to respectively minimize and maximize the value of such polynomials for a given type of input values, relying on the SB algorithm. They both return the optimal polynomial value found by the SB algorithm, along with its associated input vector.
The input types must be passed to the domain
argument:
-
spin
(default value) for a spin optimization: the optimal vector will only have ±1 values -
binary
for a binary optimization: the optimal vector will only have 0 or 1 values -
intX
for aX
-bits encoded integer optimization: the optimal vector will only have integer value encoded withX
bits or less, i.e. belonging to the range 0 to$2^{X} - 1$ .
For instance, 9-bits integer correspond to the
int9
input type and the accepted values span from 0 to 511.
import simulated_bifurcation as sb
matrix = torch.tensor([[1, 1, 2], [0, -1, -2], [-2, 0, 2]])
vector = torch.tensor([-1, 0, 2])
constant = 2.0
The package provides a polynomial
API to build quadratic multivariate polynomials from such tensors. Four options are possible.
A polynomial can be defined using coefficient tensors or SymPy expressions to define polynomials in a more natural way from mathematical equations.
The four following code snippets all create equivalent polynomials
- Using the
QuadraticPolynomial
class
from simulated_bifurcation.core import QuadraticPolynomial
With tensors
polynomial = QuadraticPolynomial(matrix, vector, constant)
With a SymPy expression
from sympy import poly, symbols
x, y, z = symbols("x y z")
expression = poly(
x**2 - y**2 + 2 * z**2
+ x * y + 2 * x * z
- 2 * y * z
- 2 * z * x
- x + 2 * z
+ 2
)
polynomial = QuadraticPolynomial(expression)
- Using the
sb.build_model
function
With tensors
polynomial = sb.build_model(matrix, vector, constant)
With a SymPy expression
from sympy import poly, symbols
x, y, z = symbols("x y z")
expression = poly(
x**2 - y**2 + 2 * z**2
+ x * y + 2 * x * z
- 2 * y * z
- 2 * z * x
- x + 2 * z
+ 2
)
polynomial = sb.build_model(expression)
The minimize
and maximize
functions allow to respectively minimize and maximize the value of such polynomials for a given type of input values, relying on the SB algorithm. They both return the optimal polynomial value found by the SB algorithm, along with its associated input vector.
# Spin minimization
spin_value, spin_vector = sb.minimize(matrix, vector, constant, domain='spin')
# Binary minimization
binary_value, binary_vector = sb.minimize(matrix, vector, constant, domain='binary')
# 3-bits integer minimization
int_value, int_vector = sb.minimize(matrix, vector, constant, domain='int3')
Or, using a SymPy expression:
# Spin minimization
spin_value, spin_vector = sb.minimize(expression, domain='spin')
# Binary minimization
binary_value, binary_vector = sb.minimize(expression, domain='binary')
# 3-bits integer minimization
int_value, int_vector = sb.minimize(expression, domain='int3')
# Spin maximization
spin_value, spin_vector = sb.maximize(matrix, vector, constant, domain='spin')
# Binary maximization
binary_value, binary_vector = sb.maximize(matrix, vector, constant, domain='binary')
# 10-bits integer maximization
int_value, int_vector = sb.maximize(matrix, vector, constant, domain='int10')
Or, using a SymPy expression:
# Spin minimization
spin_value, spin_vector = sb.maximize(expression, domain='spin')
# Binary minimization
binary_value, binary_vector = sb.maximize(expression, domain='binary')
# 3-bits integer minimization
int_value, int_vector = sb.maximize(expression, domain='int10')
For both functions, only the matrix is required, the vector and constant terms are optional.
The Simulated Bifurcation algorithm is highly parallelizable since it only relies on linear matrices equations. To take advantage of this property, this implementation offers the possibility to perform a multi-agent search of the optimal solution by evolving several spin vectors in parallel (each one being called an agent). The number of agents is set by the agents
parameter in the minimize
and maximize
functions.
💡 Tip: it is faster to run once the algorithm with N agents than to run N times the algorithm with only one agent.
# Efficient computation ✔️
sb.minimize(matrix, agents=100)
# Slower cumbersome computation ❌
for _ in range(100):
sb.minimize(matrix, agents=1)
This parallelization of the algorithm can also be utilized by performing calculations on GPUs to speed them up significantly. To do this, simply specify the calculation device
argument to cuda
when instantiating an Ising model:
sb.minimize(matrix, device='cuda')
The Simulated Bifurcation algorithm stops after a certain number of iterations, defined by the parameter max_steps
of the minimize
and maximize
functions. However, this implementation comes with the possibility to perform early stopping and save computation time by defining convergence conditions.
At regular intervals, the energy of the agents is sampled and compared with its previous value to calculate their stability period. If an agent's stability period exceeds a convergence threshold, it is considered to have converged and its value is frozen. If all agents converge before the maximum number of iterations has been reached, the algorithm stops.
- The sampling period and the convergence threshold are respectively set using the
sampling_period
andconvergence_threshold
parameters of theminimize
andmaximize
functions. - To use early stopping in the SB algorithm, set the
use_window
parameter toTrue
. - If only some agents have converged when the maximum number of iterations is reached, the algorithm stops and only these agents are considered in the results.
# Stop with maximal iterations
sb.minimize(matrix, max_steps=10000)
# Early stopping
sb.minimize(
matrix,
sampling_period=30,
convergence_threshold=50,
use_window=True,
)
By default, SB returns the best vector and objective value found. However, it is also possible to configure it so it returns all the vectors for each agent with the associated objective value. To do so, the best_only
parameter of the minimize
and maximize
functions must be set to False
(default is True
).
best_vector, best_value = sb.minimize(matrix, best_only=True)
vectors, values = sb.maximize(matrix, best_only=False)
This section deals with a more complex use of the SB algorithm, as it is closer to the quantum theory from which it is derived. To better understand the significance of the subjects at stake, we recommend reading the theory behind the SB algorithm by Goto et al..
- Goto, H., Tatsumura, K., & Dixon, A. R. (2019). Combinatorial optimization by simulating adiabatic bifurcations in nonlinear Hamiltonian systems. Science advances, 5(4), eaav2372.
- Kanao, T., & Goto, H. (2022). Simulated bifurcation assisted by thermal fluctuation. Communications Physics, 5(1), 153.
- Goto, H., Endo, K., Suzuki, M., Sakai, Y., Kanao, T., Hamakawa, Y., ... & Tatsumura, K. (2021). High-performance combinatorial optimization based on classical mechanics. Science Advances, 7(6), eabe7953.
The SB algorithm is available in four different versions (Goto et al.) that result in small variations in the algorithm general operation. The four modes are:
- Ballistic SB (bSB): uses the particles' position for the SB matrix computations; usually faster but less accurate.
- Discrete SB (dSB): uses the sign of the particles' position for the SB matrix computations; usually slower but more accurate.
- Heated ballistic SB (HbSB): uses the bSB algorithm with a supplementary non-symplectic term to allow a higher solution space exploration.
- Heated discrete SB (HdSB): uses the dSB algorithm with a supplementary non-symplectic term to allow a higher solution space exploration.
These mode can be selected setting the parameters ballistic
and heated
to either True
or False
in the Ising.optimize
method or the minimize
/maximize
functions.
sb.minimize(matrix, ballistic=True, heated=False) # bSB
sb.minimize(matrix, ballistic=False, heated=True) # HdSB
sb.maximize(matrix, ballistic=False, heated=False) # dSB
sb.maximize(matrix, ballistic=True, heated=True) # HbSB
The SB algorithm has a set of hyperparameters corresponding to physical constants derived from quantum theory, which have been fine-tuned (Goto et al.) to give the best results most of the time. Nevertheless, the relevance of specific hyperparameters may vary depending on the properties of the instances. For this purpose, the set_env
function can be used to modify their value.
# Custom hyperparameters values
sb.set_env(time_step=.1, pressure_slope=.01, heat_coefficient=.06)
# Default hyperparameters values
sb.reset_env()
A lot of mathematical problems (QUBO, Travelling Salesman Problem, MAXCUT, ...) can be written as order-two multivariate polynomials problems, and thus can be solved using the Simulated Bifurcation algorithm. Some of them are already implemented in the models
module:
🔬 Physics
- Ising model
📐 Mathematics
- Quadratic Unconstrained Binary Optimization (QUBO)
- Number partitioning
💸 Finance
- Markowitz model
You are also free to create your own models using our API. Depending on the type of model you wish to implement, you can create a subclass of the ABCModel
class to quickly and efficiently link your custom model to an Ising problem and solve it using the SB algorithm. Such a model must have a domain
class attribute that set the definition domain of all the instances.
The advantage of doing so is that your model can directly call the optimize
method that it inherits from the QuadraticPolynomial
interface without having to redefine it.
For instance, here is how the QUBO model was implemented:
The QUBO problem consists, given an upper triangular matrix
$Q$ , in finding the binary vector that minimizes the value$$\sum_{i=1}^{N} \sum_{j=1}^{N} Q_{ij}x_{i}x_{j}$$
from simulated_bifurcation.models import ABCModel
class QUBO(ABCModel):
domain = "binary"
def __init__(
self,
Q: Union[torch.Tensor, np.ndarray],
dtype: Optional[torch.dtype] = None,
device: Optional[Union[str, torch.device]] = None,
) -> None:
super().__init__(Q, dtype=dtype, device=device)
self.Q = self[2]
You can check Andrew Lucas' paper on Ising formulations of NP-complete and NP-hard problems, including all of Karp's 21 NP-complete problems.
🔎 Lucas, A. (2014). Ising formulations of many NP problems. Frontiers in physics, 2, 5.
If you are using this code for your own projects please cite our work:
@software{Porcher_Simulated_Bifurcation_SB_2023,
author = {Porcher, Francois},
month = jan,
title = {{Simulated Bifurcation (SB) algorithm for Python}},
url = {https://github.com/FrancoisPorcher/Simulated-Bifurcation},
version = {1.2.1},
year = {2024},
}