Source code for qiskit_optimization.algorithms.goemans_williamson_optimizer

# This code is part of a Qiskit project.
#
# (C) Copyright IBM 2021, 2023.
#
# This code is licensed under the Apache License, Version 2.0. You may
# obtain a copy of this license in the LICENSE.txt file in the root directory
# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
#
# Any modifications or derivative works of this code must retain this
# copyright notice, and modified files need to carry a notice indicating
# that they have been altered from the originals.

"""
Implementation of the Goemans-Williamson algorithm as an optimizer.
Requires CVXPY to run.
"""
import logging
from typing import Optional, List, Tuple, Union, cast

import numpy as np

import qiskit_optimization.optionals as _optionals
from .optimization_algorithm import (
    OptimizationResult,
    OptimizationResultStatus,
    OptimizationAlgorithm,
    SolutionSample,
)
from ..converters.flip_problem_sense import MinimizeToMaximize
from ..problems.quadratic_program import QuadraticProgram
from ..problems.variable import Variable

logger = logging.getLogger(__name__)


[docs]class GoemansWilliamsonOptimizationResult(OptimizationResult): """ Contains results of the Goemans-Williamson algorithm. The properties ``x`` and ``fval`` contain values of just one solution. Explore ``samples`` for all possible solutions. """ def __init__( self, x: Optional[Union[List[float], np.ndarray]], fval: float, variables: List[Variable], status: OptimizationResultStatus, samples: Optional[List[SolutionSample]], sdp_solution: Optional[np.ndarray] = None, ) -> None: """ Args: x: the optimal value found in the optimization. fval: the optimal function value. variables: the list of variables of the optimization problem. status: the termination status of the optimization algorithm. samples: the solution samples. sdp_solution: an SDP solution of the problem. """ super().__init__(x, fval, variables, status, samples=samples) self._sdp_solution = sdp_solution @property def sdp_solution(self) -> Optional[np.ndarray]: """ Returns: Returns an SDP solution of the problem. """ return self._sdp_solution
[docs]@_optionals.HAS_CVXPY.require_in_instance class GoemansWilliamsonOptimizer(OptimizationAlgorithm): """ Goemans-Williamson algorithm to approximate the max-cut of a problem. The quadratic program for max-cut is given by: max sum_{i,j<i} w[i,j]*x[i]*(1-x[j]) Therefore the quadratic term encodes the negative of the adjacency matrix of the graph. """ def __init__( self, num_cuts: int, sort_cuts: bool = True, unique_cuts: bool = True, seed: int = 0, ): """ Args: num_cuts: Number of cuts to generate. sort_cuts: True if sort cuts by their values. unique_cuts: The solve method returns only unique cuts, thus there may be less cuts than ``num_cuts``. seed: A seed value for the random number generator. """ super().__init__() self._num_cuts = num_cuts self._sort_cuts = sort_cuts self._unique_cuts = unique_cuts np.random.seed(seed)
[docs] def get_compatibility_msg(self, problem: QuadraticProgram) -> str: """Checks whether a given problem can be solved with the optimizer implementing this method. Args: problem: The optimization problem to check compatibility. Returns: Returns the incompatibility message. If the message is empty no issues were found. """ message = "" if problem.get_num_binary_vars() != problem.get_num_vars(): message = ( f"Only binary variables are supported, while the total number of variables " f"{problem.get_num_vars()} and there are {problem.get_num_binary_vars()} " f"binary variables across them" ) return message
[docs] def solve(self, problem: QuadraticProgram) -> OptimizationResult: """ Returns a list of cuts generated according to the Goemans-Williamson algorithm. Args: problem: The quadratic problem that encodes the max-cut problem. Returns: cuts: A list of generated cuts. """ # pylint: disable=import-error from cvxpy import DCPError, DGPError, SolverError self._verify_compatibility(problem) min2max = MinimizeToMaximize() problem = min2max.convert(problem) adj_matrix = self._extract_adjacency_matrix(problem) try: chi = self._solve_max_cut_sdp(adj_matrix) except (DCPError, DGPError, SolverError): logger.error("Can't solve SDP problem") return GoemansWilliamsonOptimizationResult( x=[], fval=0, variables=problem.variables, status=OptimizationResultStatus.FAILURE, samples=[], ) cuts = self._generate_random_cuts(chi, len(adj_matrix)) numeric_solutions = [ (cuts[i, :], self.max_cut_value(cuts[i, :], adj_matrix)) for i in range(self._num_cuts) ] if self._sort_cuts: numeric_solutions.sort(key=lambda x: -x[1]) if self._unique_cuts: numeric_solutions = self._get_unique_cuts(numeric_solutions) numeric_solutions = numeric_solutions[: self._num_cuts] samples = [ SolutionSample( x=solution[0], fval=solution[1], probability=1.0 / len(numeric_solutions), status=OptimizationResultStatus.SUCCESS, ) for solution in numeric_solutions ] return cast( GoemansWilliamsonOptimizationResult, self._interpret( x=samples[0].x, problem=problem, converters=[min2max], result_class=GoemansWilliamsonOptimizationResult, samples=samples, ), )
def _get_unique_cuts( self, solutions: List[Tuple[np.ndarray, float]] ) -> List[Tuple[np.ndarray, float]]: """ Returns: Unique Goemans-Williamson cuts. """ # Remove symmetry in the cuts to chose the unique ones. # Cuts 010 and 101 are symmetric(same cut), so we convert all cuts # starting from 1 to start from 0. In the next loop repetitive cuts will be removed. for idx, cut in enumerate(solutions): if cut[0][0] == 1: solutions[idx] = ( np.array([0 if _ == 1 else 1 for _ in cut[0]]), cut[1], ) seen_cuts = set() unique_cuts = [] for cut in solutions: cut_str = "".join([str(_) for _ in cut[0]]) if cut_str in seen_cuts: continue seen_cuts.add(cut_str) unique_cuts.append(cut) return unique_cuts @staticmethod def _extract_adjacency_matrix(problem: QuadraticProgram) -> np.ndarray: """ Extracts the adjacency matrix from the given quadratic program. Args: problem: A QuadraticProgram describing the max-cut optimization problem. Returns: adjacency matrix of the graph. """ adj_matrix = -problem.objective.quadratic.coefficients.toarray() adj_matrix = (adj_matrix + adj_matrix.T) / 2 return adj_matrix def _solve_max_cut_sdp(self, adj_matrix: np.ndarray) -> np.ndarray: """ Calculates the maximum weight cut by generating |V| vectors with a vector program, then generating a random plane that cuts the vertices. This is the Goemans-Williamson algorithm that gives a .878-approximation. Returns: chi: a list of length |V| where the i-th element is +1 or -1, representing which set the it-h vertex is in. Returns None if an error occurs. """ # pylint: disable=import-error import cvxpy as cvx num_vertices = len(adj_matrix) constraints, expr = [], 0 # variables x = cvx.Variable((num_vertices, num_vertices), PSD=True) # constraints for i in range(num_vertices): constraints.append(x[i, i] == 1) # objective function expr = cvx.sum(cvx.multiply(adj_matrix, (np.ones((num_vertices, num_vertices)) - x))) # solve problem = cvx.Problem(cvx.Maximize(expr), constraints) problem.solve() return x.value def _generate_random_cuts(self, chi: np.ndarray, num_vertices: int) -> np.ndarray: """ Random hyperplane partitions vertices. Args: chi: a list of length |V| where the i-th element is +1 or -1, representing which set the i-th vertex is in. num_vertices: the number of vertices in the graph Returns: An array of random cuts. """ eigenvalues = np.linalg.eigh(chi)[0] if min(eigenvalues) < 0: chi = chi + (1.001 * abs(min(eigenvalues)) * np.identity(num_vertices)) elif min(eigenvalues) == 0: chi = chi + 0.00001 * np.identity(num_vertices) x = np.linalg.cholesky(chi).T r = np.random.normal(size=(self._num_cuts, num_vertices)) return (np.dot(r, x) > 0) + 0
[docs] @staticmethod def max_cut_value(x: np.ndarray, adj_matrix: np.ndarray): """Compute the value of a cut from an adjacency matrix and a list of binary values. Args: x: a list of binary value in numpy array. adj_matrix: adjacency matrix. Returns: float: value of the cut. """ cut_matrix = np.outer(x, (1 - x)) return np.sum(adj_matrix * cut_matrix)