Use correct local cost to decide on switching

src/lava/lib/optimization/problems/problems.py

@@ -65,8 +65,8 @@ def __init__(self, q: npt.ArrayLike):
 
         Parameters
         ----------
-        q: squared Q matrix defining the QUBO problem over a binary
-        vector x as: minimize x^T*Q*x.
+        q: squared, symmetric, int Q matrix defining the QUBO problem over a
+        binary vector x as: minimize x^T*Q*x.
         """
         super().__init__()
         self.validate_input(q)
@@ -109,20 +109,26 @@ def evaluate_cost(self, solution: np.ndarray) -> int:
         return int(self._q_cost(solution))
 
     def validate_input(self, q):
-        """Validate the cost coefficient is a square matrix.
+        """Validate that cost coefficient is a square, symmetric, int matrix.
 
         Parameters
         ----------
         q: Quadratic coefficient of the cost function.
 
         """
+
         m, n = q.shape
         if m != n:
             raise ValueError("q matrix is not a square matrix.")
         if not issubclass(q.dtype.type, np.integer):
             raise NotImplementedError(
                 "Non integer q matrices are not supported yet."
             )
+        # matrix must be symmetric for current implementation
+        if not np.allclose(q, q.T, rtol=1e-05, atol=1e-08):
+            raise NotImplementedError(
+                "Only symmetric matrixes are currently supported."
+            )
 
     def verify_solution(self, solution):
         raise NotImplementedError

src/lava/lib/optimization/solvers/generic/hierarchical_processes.py

@@ -9,7 +9,6 @@
 from lava.magma.core.process.ports.ports import InPort, OutPort
 from lava.magma.core.process.process import AbstractProcess, LogConfig
 from lava.magma.core.process.variable import Var
-from numpy import typing as npt
 
 
 class ContinuousVariablesProcess(AbstractProcess):
@@ -21,7 +20,7 @@ def __init__(
         shape: ty.Tuple[int, ...],
         problem: OptimizationProblem,
         backend,
-        hyperparameters: ty.Dict[str, ty.Union[int, npt.ArrayLike]] = None,
+        hyperparameters: ty.Dict[str, ty.Union[int, npty.ArrayLike]] = None,
         name: ty.Optional[str] = None,
         log_config: ty.Optional[LogConfig] = None,
     ) -> None:
@@ -67,8 +66,8 @@ class DiscreteVariablesProcess(AbstractProcess):
     def __init__(
         self,
         shape: ty.Tuple[int, ...],
-        cost_diagonal: npt.ArrayLike = None,
-        hyperparameters: ty.Dict[str, ty.Union[int, npt.ArrayLike]] = None,
+        cost_diagonal: npty.ArrayLike = None,
+        hyperparameters: ty.Dict[str, ty.Union[int, npty.ArrayLike]] = None,
         name: ty.Optional[str] = None,
         log_config: ty.Optional[LogConfig] = None,
     ) -> None:
@@ -77,7 +76,7 @@ def __init__(
         ----------
         shape: tuple
             A tuple of the form (number of variables, domain size).
-        cost_diagonal: npt.ArrayLike
+        cost_diagonal: npty.ArrayLike
             The diagonal of the coefficient of the quadratic term on the cost
             function.
         hyperparameters: dict, optional
@@ -180,7 +179,7 @@ def __init__(
         shape_out: ty.Tuple[int, ...],
         problem: OptimizationProblem,
         backend,
-        hyperparameters: ty.Dict[str, ty.Union[int, npt.ArrayLike]] = None,
+        hyperparameters: ty.Dict[str, ty.Union[int, npty.ArrayLike]] = None,
         name: ty.Optional[str] = None,
         log_config: ty.Optional[LogConfig] = None,
     ) -> None:
@@ -246,36 +245,36 @@ class StochasticIntegrateAndFire(AbstractProcess):
     def __init__(
         self,
         *,
-        step_size: npt.ArrayLike,
+        step_size: npty.ArrayLike,
         shape: ty.Tuple[int, ...] = (1,),
-        init_state: npt.ArrayLike = 0,
-        noise_amplitude: npt.ArrayLike = 1,
-        noise_precision: npt.ArrayLike = 8,
-        sustained_on_tau: npt.ArrayLike = -3,
-        threshold: npt.ArrayLike = 10,
-        cost_diagonal: npt.ArrayLike = 0,
+        init_state: npty.ArrayLike = 0,
+        noise_amplitude: npty.ArrayLike = 1,
+        noise_precision: npty.ArrayLike = 8,
+        sustained_on_tau: npty.ArrayLike = -3,
+        threshold: npty.ArrayLike = 10,
+        cost_diagonal: npty.ArrayLike = 0,
         name: ty.Optional[str] = None,
         log_config: ty.Optional[LogConfig] = None,
-        init_value: npt.ArrayLike = 0,
+        init_value: npty.ArrayLike = 0,
     ) -> None:
         """
         Parameters
         ----------
         shape: tuple
             The shape of the set of dynamical systems to be created.
-        init_state: npt.ArrayLike, optional
+        init_state: npty.ArrayLike, optional
             The starting value of the state variable.
-        step_size: npt.ArrayLike, optional
+        step_size: npty.ArrayLike, optional
             a value to be added to the state variable at each timestep.
-        noise_amplitude: npt.ArrayLike, optional
+        noise_amplitude: npty.ArrayLike, optional
             The width/range for the stochastic perturbation to the state
             variable. A random number within this range will be added to the
             state variable at each timestep.
-        steps_to_fire: npt.ArrayLike, optional
+        steps_to_fire: npty.ArrayLike, optional
             After how many timesteps would the dynamical system fire and reset
             without stochastic perturbation. Note that if noise_amplitude > 0,
             the system will stochastically deviate from this value.
-        cost_diagonal: npt.ArrayLike, optional
+        cost_diagonal: npty.ArrayLike, optional
             The linear coefficients on the cost function of the optimization
             problem where this system will be used.
         name: str, optional
@@ -342,45 +341,45 @@ class NEBMAbstract(AbstractProcess):
     def __init__(
         self,
         *,
-        temperature: npt.ArrayLike,
+        temperature: npty.ArrayLike,
         refract: ty.Optional[ty.Union[int, npty.NDArray]],
         refract_counter: ty.Optional[ty.Union[int, npty.NDArray]],
         shape: ty.Tuple[int, ...] = (1,),
-        init_state: npt.ArrayLike = 0,
-        input_duration: npt.ArrayLike = 6,
-        min_state: npt.ArrayLike = 1000,
-        min_integration: npt.ArrayLike = -1000,
-        cost_diagonal: npt.ArrayLike = 0,
+        init_state: npty.ArrayLike = 0,
+        input_duration: npty.ArrayLike = 6,
+        min_state: npty.ArrayLike = 1000,
+        min_integration: npty.ArrayLike = -1000,
+        cost_diagonal: npty.ArrayLike = 0,
         name: ty.Optional[str] = None,
         log_config: ty.Optional[LogConfig] = None,
-        init_value: npt.ArrayLike = 0,
+        init_value: npty.ArrayLike = 0,
     ) -> None:
         """
 
         Parameters
         ----------
         shape: tuple
             The shape of the set of dynamical systems to be created.
-        init_state: npt.ArrayLike, optional
+        init_state: npty.ArrayLike, optional
             The starting value of the state variable.
-        temperature: npt.ArrayLike, optional
+        temperature: npty.ArrayLike, optional
             the temperature of the systems, defining the level of noise.
-        input_duration: npt.ArrayLike, optional
+        input_duration: npty.ArrayLike, optional
             Number of timesteps by which each input should be preserved.
-        min_state: npt.ArrayLike, optional
+        min_state: npty.ArrayLike, optional
             The minimum value for the state variable. The state variable will be
             truncated at this value if updating results in a lower value.
-        min_integration: npt.ArrayLike, optional
+        min_integration: npty.ArrayLike, optional
             The minimum value for the total input (addition of all valid inputs
             at a given timestep). The total input value will be truncated at
             this value if adding current and preserved inputs results in a lower
             value.
-        refract: npt.ArrayLike, optional
+        refract: npty.ArrayLike, optional
             Number of timesteps to wait after firing and reset before resuming
             updating.
-        refract_counter: npt.ArrayLike, optional
+        refract_counter: npty.ArrayLike, optional
             Number of timesteps to initially suppress a unit firing.
-        cost_diagonal: npt.ArrayLike, optional
+        cost_diagonal: npty.ArrayLike, optional
             The linear coefficients on the cost function of the optimization
             problem where this system will be used.
         name: str, optional
@@ -450,6 +449,7 @@ class NEBMSimulatedAnnealingAbstract(AbstractProcess):
     def __init__(
         self,
         *,
+        cost_diagonal: npty.ArrayLike,
         max_temperature: int = 10,
         min_temperature: int = 0,
         delta_temperature: int = 1,
@@ -458,12 +458,11 @@ def __init__(
         refract_scaling: int = 14,
         refract: ty.Optional[ty.Union[int, npty.NDArray]],
         shape: ty.Tuple[int, ...] = (1,),
-        init_state: npt.ArrayLike = 0,
-        min_integration: npt.ArrayLike = -1000,
-        cost_diagonal: npt.ArrayLike = 0,
+        init_state: npty.ArrayLike = 0,
+        min_integration: npty.ArrayLike = -1000,
         name: ty.Optional[str] = None,
         log_config: ty.Optional[LogConfig] = None,
-        init_value: npt.ArrayLike = 0,
+        init_value: npty.ArrayLike = 0,
         annealing_schedule: str = 'linear',
         neuron_model: str,
     ) -> None:
@@ -473,19 +472,19 @@ def __init__(
         ----------
         shape: tuple
             The shape of the set of dynamical systems to be created.
-        init_state: npt.ArrayLike, optional
+        init_state: npty.ArrayLike, optional
             The starting value of the state variable.
-        temperature: npt.ArrayLike, optional
+        temperature: npty.ArrayLike, optional
             the temperature of the systems, defining the level of noise.
-        min_integration: npt.ArrayLike, optional
+        min_integration: npty.ArrayLike, optional
             The minimum value for the total input (addition of all valid inputs
             at a given timestep). The total input value will be truncated at
             this value if adding current and preserved inputs results in a lower
             value.
-        refractory_period: npt.ArrayLike, optional
+        refractory_period: npty.ArrayLike, optional
             Number of timesteps to wait after firing and reset before resuming
             updating.
-        cost_diagonal: npt.ArrayLike, optional
+        cost_diagonal: npty.ArrayLike, optional
             The linear coefficients on the cost function of the optimization
             problem where this system will be used.
         name: str, optional

src/lava/lib/optimization/solvers/generic/nebm/process.py

@@ -70,6 +70,7 @@ class NEBMSimulatedAnnealing(AbstractProcess):
     def __init__(
         self,
         *,
+        cost_diagonal: npty.ArrayLike,
         shape: ty.Tuple[int, ...],
         max_temperature: int,
         min_temperature: int,
@@ -109,6 +110,7 @@ def __init__(
             delta_temperature=delta_temperature,
             steps_per_temperature=steps_per_temperature,
             refract=refract,
+            cost_diagonal=cost_diagonal,
             refract_scaling=refract_scaling,
             exp_temperature=exp_temperature,
             neuron_model=neuron_model,

src/lava/lib/optimization/solvers/generic/solution_finder/models.py

@@ -192,4 +192,9 @@ def _get_init_state(
         else:
             q_diag = cost_coefficients[2].init.diagonal()
 
-        return q_off_diag @ init_value + q_diag
+        # nebm-sa-refract neuron models have a different definition of state
+        neuron_model = hyperparameters.get("neuron_model", "nebm")
+        if neuron_model == 'nebm-sa-refract':
+            return q_off_diag @ init_value
+        else:
+            return q_off_diag @ init_value + q_diag

src/lava/lib/optimization/solvers/generic/solver.py

@@ -6,6 +6,7 @@
 import numpy as np
 
 from dataclasses import dataclass
+import warnings
 
 from lava.lib.optimization.solvers.generic.qp.models import (
     PyPIneurPIPGeqModel,
@@ -33,6 +34,7 @@
 from lava.proc.monitor.process import Monitor
 from lava.utils.profiler import Profiler
 
+from lava.lib.optimization.problems.problems import QUBO
 from lava.lib.optimization.problems.problems import OptimizationProblem
 from lava.lib.optimization.solvers.generic.builder import SolverProcessBuilder
 from lava.lib.optimization.solvers.generic.cost_integrator.process import (
@@ -271,6 +273,7 @@ def solve(self, config: SolverConfig = SolverConfig()) -> SolverReport:
         report: SolverReport
             An object containing all the data generated by the execution.
         """
+        self._validate_input_solve(config)
         run_condition, run_cfg = self._prepare_solver(config)
         self.solver_process.run(condition=run_condition, run_cfg=run_cfg)
         best_state, best_cost, best_timestep = self._get_results(config)
@@ -535,3 +538,26 @@ def _get_and_decode_continuous_vars(self, idx: int):
             self.solver_process.finders[idx].variables_assignment.get()
         )
         return solution
+
+    def _validate_input_solve(self, config: SolverConfig) -> None:
+        """
+        Validates the user input provided to the solve() method of the
+        OptimizationSolver.
+
+        Parameters
+        ----------
+        config: SolverConfig
+            Solver configuration used.
+
+        Raises
+        ----------
+        """
+
+        if (config.backend in CPUS) and isinstance(self.problem, QUBO):
+            warnings.warn(f"Please note that we are currently advancing only "
+                          f"the Loihi 2 backend of the QUBO solver. Until "
+                          f"further notice, the CPU implementation uses "
+                          f"an outdated algorithm. We thus recommend running "
+                          f"your workloads on Loihi 2 to leverage the latest "
+                          f"algorithms for optimal latency, energy "
+                          f"consumption, and solution accuracy. ")

src/lava/lib/optimization/solvers/generic/sub_process_models.py

@@ -314,6 +314,7 @@ def __init__(self, proc):
                 steps_per_temperature=steps_per_temperature,
                 refract=refract,
                 refract_scaling=refract_scaling,
+                cost_diagonal=diagonal,
                 init_value=init_value,
                 init_state=init_state,
                 neuron_model=neuron_model,
@@ -457,6 +458,7 @@ class NEBMSimulatedAnnealingAbstractModel(AbstractSubProcessModel):
 
     def __init__(self, proc):
         shape = proc.proc_params.get("shape", (1,))
+        cost_diagonal = proc.proc_params.get("cost_diagonal")
         max_temperature = proc.proc_params.get("max_temperature", 10)
         min_temperature = proc.proc_params.get("min_temperature", 0)
         delta_temperature = proc.proc_params.get("delta_temperature", 1)
@@ -477,6 +479,7 @@ def __init__(self, proc):
             exp_temperature=exp_temperature,
             steps_per_temperature=steps_per_temperature,
             refract_scaling=refract_scaling,
+            cost_diagonal=cost_diagonal,
             refract=refract,
             init_value=init_value,
             init_state=init_state,

tests/lava/lib/optimization/utils/test_report_analyzer.py

@@ -15,9 +15,9 @@ def prepare_problem_and_report():
     np.random.seed(0)
     size = 5
     timeout = 100
-    problem = QUBO(
-        q=np.random.randint(0, 20, size=(size, size), dtype=np.int32)
-    )
+    q = np.random.randint(0, 20, size=(size, size), dtype=np.int32)
+    q_symm = ((q + q.T) / 2).astype(int)
+    problem = QUBO(q=q_symm)
     states = np.random.randint(0, 2, size=(timeout, size), dtype=np.int32)
     costs = list(map(problem.evaluate_cost, states))
     report = SolverReport(