Source code for qiskit_optimization.algorithms.qrao.quantum_random_access_encoding

# This code is part of a Qiskit project.
#
# (C) Copyright IBM 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.

"""The Quantum Random Access Encoding module."""
from __future__ import annotations

from collections import defaultdict
from functools import reduce
from typing import cast

import networkx as nx
import numpy as np
from qiskit import QuantumCircuit
from qiskit.quantum_info import SparsePauliOp

from qiskit_optimization.exceptions import QiskitOptimizationError
from qiskit_optimization.problems.quadratic_program import QuadraticProgram


def _z_to_31p_qrac_basis_circuit(bases: list[int], bit_flip: int = 0) -> QuantumCircuit:
    """Return the circuit that implements the rotation to the (3,1,p)-QRAC.

    Args:
        bases: The basis, 0, 1, 2, or 3, for the qubit.
        bit_flip: Whether to flip the state of the qubit. 1 for flip, 0 for no flip.

    Returns:
        The ``QuantumCircuit`` implementing the rotation to the (3,1,p)-QRAC.

    Raises:
        ValueError: If the basis is not 0, 1, 2, or 3
    """
    circ = QuantumCircuit(len(bases))
    BETA = np.arccos(1 / np.sqrt(3))  # pylint: disable=invalid-name

    for i, base in enumerate(reversed(bases)):
        if bit_flip:
            # if bit_flip == 1: then flip the state of the qubit to |1>
            circ.x(i)

        if base == 0:
            circ.r(-BETA, -np.pi / 4, i)
        elif base == 1:
            circ.r(np.pi - BETA, np.pi / 4, i)
        elif base == 2:
            circ.r(np.pi + BETA, np.pi / 4, i)
        elif base == 3:
            circ.r(BETA, -np.pi / 4, i)
        else:
            raise ValueError(f"Unknown basis: {base}. Basis must be 0, 1, 2, or 3.")
    return circ


def _z_to_21p_qrac_basis_circuit(bases: list[int], bit_flip: int = 0) -> QuantumCircuit:
    """Return the circuit that implements the rotation to the (2,1,p)-QRAC.

    Args:
        bases: The basis, 0, 1, for the qubit.
        bit_flip: Whether to flip the state of the qubit. 1 for flip, 0 for no flip.

    Returns:
        The ``QuantumCircuit`` implementing the rotation to the (2,1,p)-QRAC.

    Raises:
        ValueError: if the basis is not 0 or 1
    """
    circ = QuantumCircuit(len(bases))

    for i, base in enumerate(reversed(bases)):
        if bit_flip:
            # if bit_flip == 1: then flip the state of the qubit to |1>
            circ.x(i)

        if base == 0:
            circ.r(-1 * np.pi / 4, -np.pi / 2, i)
        elif base == 1:
            circ.r(-3 * np.pi / 4, -np.pi / 2, i)
        else:
            raise ValueError(f"Unknown basis: {base}. Basis must be 0, 1.")
    return circ


def _qrac_state_prep_1q(bit_list: list[int]) -> QuantumCircuit:
    """
    Return the circuit that prepares the state for a (1,1,p), (2,1,p), or (3,1,p)-QRAC.

    Args:
        bit_list: The bitstring to prepare. If 1 argument is given, then a (1,1,p)-QRAC is generated.
            If 2 arguments are given, then a (2,1,p)-QRAC is generated. If 3 arguments are given,
            then a (3,1,p)-QRAC is generated.

    Returns:
        The ``QuantumCircuit`` implementing the state preparation.

    Raises:
        TypeError: if the number of arguments is not 1, 2, or 3
        ValueError: if any of the arguments are not 0 or 1
    """
    # pylint: disable=C0401
    if len(bit_list) not in (1, 2, 3):
        raise TypeError(f"qrac_state_prep_1q requires 1, 2, or 3 arguments, not {len(bit_list)}.")
    if not all(bit in (0, 1) for bit in bit_list):
        raise ValueError("Each argument to qrac_state_prep_1q must be 0 or 1.")

    if len(bit_list) == 3:
        # Prepare (3,1,p)-qrac
        # In the following lines, the input bits are XORed to match the
        # conventions used in the paper.
        # To understand why this transformation happens,
        # observe that the two states that define each magic basis
        # correspond to the same bitstrings but with a global bit flip.
        # Thus the three bits of information we use to construct these states are:
        # base_index0,base_index1 : two bits to pick one of four magic bases
        # bit_flip: one bit to indicate which magic basis projector we are interested in.

        bit_flip = bit_list[0] ^ bit_list[1] ^ bit_list[2]
        base_index0 = bit_list[1] ^ bit_list[2]
        base_index1 = bit_list[0] ^ bit_list[2]

        # This is a convention chosen to be consistent with https://arxiv.org/abs/2111.03167
        # See SI:4 second paragraph and observe that π+ = |0X0|, π- = |1X1|
        base = [2 * base_index0 + base_index1]
        circ = _z_to_31p_qrac_basis_circuit(base, bit_flip)

    elif len(bit_list) == 2:
        # Prepare (2,1,p)-qrac
        # (00,01) or (10,11)
        bit_flip = bit_list[0]
        # (00,11) or (01,10)
        base_index0 = bit_list[0] ^ bit_list[1]

        # This is a convention chosen to be consistent with https://arxiv.org/abs/2111.03167
        # # See SI:4 second paragraph and observe that π+ = |0X0|, π- = |1X1|
        base = [base_index0]
        circ = _z_to_21p_qrac_basis_circuit(base, bit_flip)

    else:
        bit_flip = bit_list[0]
        circ = QuantumCircuit(1)
        if bit_flip:
            circ.x(0)

    return circ


def _qrac_state_prep_multi_qubit(
    x: list[int],
    q2vars: list[list[int]],
    max_vars_per_qubit: int,
) -> QuantumCircuit:
    """Prepares a multi qubit QRAC state.

    Args:
        x: The state of each decision variable (0 or 1).
        q2vars: A list of lists of integers. Each inner list contains the indices of decision variables
        mapped to a specific qubit.
        max_vars_per_qubit: The maximum number of decision variables that can be mapped to a
        single qubit.

    Returns:
        A QuantumCircuit object representing the prepared state.

    Raises:
        ValueError: If any qubit is associated with more than `max_vars_per_qubit` variables.
        ValueError: If a decision variable in ``q2vars`` is not included in `x`.
        ValueError: If there are unused decision variables in `x` after mapping to qubits.
    """
    # Create a set of all remaining decision variables
    remaining_dvars = set(range(len(x)))
    # Create a list to store the binary mappings of each qubit to its corresponding decision variables
    variable_mappings: list[list[int]] = []
    # Check that each qubit is associated with at most max_vars_per_qubit variables
    for qi_vars in q2vars:
        if len(qi_vars) > max_vars_per_qubit:
            raise ValueError(
                "Each qubit is expected to be associated with at most "
                f"`max_vars_per_qubit` ({max_vars_per_qubit}) variables, "
                f"not {len(qi_vars)} variables."
            )
        # Create a list to store the binary mapping of the current qubit
        qi_bits: list[int] = []

        # Map each decision variable associated with the current qubit to a binary value and add it
        # to the qubit bits
        for dvar in qi_vars:
            try:
                qi_bits.append(x[dvar])
            except IndexError:
                raise ValueError(f"Decision variable not included in dvars: {dvar}") from None
            try:
                remaining_dvars.remove(dvar)
            except KeyError:
                raise ValueError(
                    f"Unused decision variable(s) in dvars: {remaining_dvars}"
                ) from None

        # Pad with zeros if necessary
        while len(qi_bits) < max_vars_per_qubit:
            qi_bits.append(0)

        variable_mappings.append(qi_bits)

    # Raise an error if not all decision variables are used
    if remaining_dvars:
        raise ValueError(f"Not all dvars were included in q2vars: {remaining_dvars}")

    # Prepare the individual qrac circuit and combine them into a multi qubit circuit
    qracs = [_qrac_state_prep_1q(qi_bits) for qi_bits in variable_mappings]
    qrac_circ = reduce(lambda x, y: x.tensor(y), qracs)
    return qrac_circ


[docs]class QuantumRandomAccessEncoding: """This class specifies a Quantum Random Access Code that can be used to encode the binary variables of a QUBO (quadratic unconstrained binary optimization problem). """ # This defines the convention of the Pauli operators (and their ordering) # for each encoding. _OPERATORS = ( (SparsePauliOp("Z"),), # (1,1,1) QRAC (SparsePauliOp("X"), SparsePauliOp("Z")), # (2,1,p) QRAC, p ≈ 0.85 (SparsePauliOp("X"), SparsePauliOp("Y"), SparsePauliOp("Z")), # (3,1,p) QRAC, p ≈ 0.79 ) def __init__(self, max_vars_per_qubit: int = 3): """ Args: max_vars_per_qubit: The maximum number of decision variables per qubit. Integer values 1, 2 and 3 are supported (default to 3). """ if max_vars_per_qubit not in (1, 2, 3): raise ValueError("max_vars_per_qubit must be 1, 2, or 3") self._ops = self._OPERATORS[max_vars_per_qubit - 1] self._qubit_op: SparsePauliOp | None = None self._offset: float | None = None self._problem: QuadraticProgram | None = None self._var2op: dict[int, tuple[int, SparsePauliOp]] = {} self._q2vars: list[list[int]] = [] self._frozen = False @property def num_qubits(self) -> int: """Number of qubits""" return len(self._q2vars) @property def num_vars(self) -> int: """Number of decision variables""" return len(self._var2op) @property def max_vars_per_qubit(self) -> int: """Maximum number of variables per qubit""" return len(self._ops) @property def var2op(self) -> dict[int, tuple[int, SparsePauliOp]]: """Maps each decision variable to ``(qubit_index, operator)``""" return self._var2op @property def q2vars(self) -> list[list[int]]: """Each element contains the list of decision variable indices encoded on that qubit""" return self._q2vars @property def compression_ratio(self) -> float: """Compression ratio. Number of decision variables divided by number of qubits""" return self.num_vars / self.num_qubits @property def minimum_recovery_probability(self) -> float: """Minimum recovery probability, as set by ``max_vars_per_qubit``""" n = self.max_vars_per_qubit return (1 + 1 / np.sqrt(n)) / 2 @property def qubit_op(self) -> SparsePauliOp: """Relaxed Hamiltonian operator. Raises: RuntimeError: If the objective function has not been set yet. Use the ``encode`` method to construct the Hamiltonian, or make sure that the objective function has been set. """ if self._qubit_op is None: raise RuntimeError( "Cannot return the relaxed Hamiltonian operator: no objective function has been " "provided yet. Use the ``encode`` method to construct the Hamiltonian, or make " "sure that the objective function has been set." ) return self._qubit_op @property def offset(self) -> float: """Relaxed Hamiltonian offset Raises: RuntimeError: If the offset has not been set yet. Use the ``encode`` method to construct the Hamiltonian, or make sure that the objective function has been set. """ if self._offset is None: raise RuntimeError( "Cannot return the relaxed Hamiltonian offset: The offset attribute cannot be " "accessed until the ``encode`` method has been called to generate the qubit " "Hamiltonian. Please call ``encode`` first." ) return self._offset @property def problem(self) -> QuadraticProgram: """The ``QuadraticProgram`` encoding a QUBO optimization problem Raises: RuntimeError: If the ``QuadraticProgram`` has not been set yet. Use the ``encode`` method to set the problem. """ if self._problem is None: raise RuntimeError( "This object has not been associated with a ``QuadraticProgram``. " "Please use the ``encode`` method to set the problem." ) return self._problem
[docs] def freeze(self): """Freeze the object to prevent further modification. Once an instance of this class is frozen, ``encode`` can no longer be called. """ if not self._frozen: self._qubit_op = self._qubit_op.simplify(atol=0) self._frozen = True
@property def frozen(self) -> bool: """Whether the object is frozen or not.""" return self._frozen def _add_variables(self, variables: list[int]) -> None: """Add variables to the Encoding object. Args: variables: A list of variable indices to be added. Raises: ValueError: If added variables are not unique. ValueError: If added variables collide with existing ones. """ # NOTE: If this is called multiple times, it *always* adds an # additional qubit (see final line), even if aggregating them into a # single call would have resulted in fewer qubits. # Check if variables is empty if not variables: return # Check if variables are unique if len(variables) != len(set(variables)): raise ValueError("Added variables must be unique") # Check if variables collide with existing ones for v in variables: if v in self._var2op: raise ValueError("Added variables cannot collide with existing ones") # Calculate the number of new qubits required for the added variables. n = len(self._ops) old_num_qubits = len(self._q2vars) num_new_qubits = int(np.ceil(len(variables) / n)) # Add the new qubits to _q2vars. for _ in range(num_new_qubits): self._q2vars.append([]) # Associate each added variable with a qubit and operator. for i, v in enumerate(variables): qubit, op = divmod(i, n) qubit_index = old_num_qubits + qubit self._var2op[v] = (qubit_index, self._ops[op]) self._q2vars[qubit_index].append(v) def _add_term(self, w: float, *variables: int) -> None: """Add a term to the Hamiltonian. Args: weight: the coefficient for the term *variables: the list of variables for the term """ # Eq. (31) in https://arxiv.org/abs/2111.03167 assumes a weight-2 # Pauli operator. To generalize, we replace the `d` in that equation # with `d_prime`, defined as follows: d_prime = np.sqrt(self.max_vars_per_qubit) ** len(variables) op = self._term2op(*variables) * (w * d_prime) if w == 0.0: return if self._qubit_op is None: self._qubit_op = op else: self._qubit_op += op def _term2op(self, *variables: int) -> SparsePauliOp: """Construct a ``SparsePauliOp`` that is a tensor product of encoded decision variables. Args: *variables: The indices of the decision variables to encode. Returns: The encoded ``SparsePauliOp`` representing the product of the provided variables. Raises: QiskitOptimizationError: If any of the decision variables to be encoded collide in qubit space. """ ops = [SparsePauliOp("I")] * self.num_qubits done = set() for x in variables: pos, op = self._var2op[x] if pos in done: raise QiskitOptimizationError(f"Collision of variables: {variables}") ops[pos] = op done.add(pos) pauli_op = reduce(lambda x, y: x.tensor(y), ops) return pauli_op @staticmethod def _generate_ising_coefficients( problem: QuadraticProgram, ) -> tuple[float, np.ndarray, np.ndarray]: """Generate coefficients of Hamiltonian from a given problem.""" num_vars = problem.get_num_vars() # set a sign corresponding to a maximized or minimized problem: # 1 is for minimized problem, -1 is for maximized problem. sense = problem.objective.sense.value # convert a constant part of the objective function into Hamiltonian. offset = problem.objective.constant * sense # convert linear parts of the objective function into Hamiltonian. linear = np.zeros(num_vars) for idx, coef in problem.objective.linear.to_dict().items(): idx = cast(int, idx) weight = coef * sense / 2 linear[idx] -= weight offset += weight # convert quadratic parts of the objective function into Hamiltonian. quad = np.zeros((num_vars, num_vars)) for (i, j), coef in problem.objective.quadratic.to_dict().items(): i = cast(int, i) j = cast(int, j) weight = coef * sense / 4 if i == j: linear[i] -= 2 * weight offset += 2 * weight else: quad[i, j] += weight linear[i] -= weight linear[j] -= weight offset += weight return offset, linear, quad @staticmethod def _find_variable_partition(quad: np.ndarray) -> dict[int, list[int]]: """Find the variable partition of the quad based on the node coloring of the graph Args: coefficients of the quadratic part of the Hamiltonian. Returns: A dictionary of the variable partition of the quad based on the node coloring. """ # pylint: disable=E1101 color2node: dict[int, list[int]] = defaultdict(list) num_nodes = quad.shape[0] graph = nx.Graph() graph.add_nodes_from(range(num_nodes)) graph.add_edges_from(list(zip(*np.where(quad != 0)))) node2color = nx.greedy_color(graph) for node, color in sorted(node2color.items()): color2node[color].append(node) return color2node
[docs] def encode(self, problem: QuadraticProgram) -> None: """ Encodes a given ``QuadraticProgram`` as a (n,1,p) Quantum Random Access Code (QRAC) relaxed Hamiltonian. It accomplishes this by mapping each binary decision variable to one qubit of the QRAC. The encoding is designed to ensure that the problem's objective function commutes with the QRAC encoding. After the function is called, it sets the following attributes: - qubit_op: The qubit operator that encodes the input ``QuadraticProgram``. - offset: The constant value in the encoded Hamiltonian. - problem: The original ``QuadraticProgram`` used for encoding. Inputs: problem: A ``QuadraticProgram`` encoding a QUBO optimization problem Raises: QiskitOptimizationError: If this method is called more than once on the same object. QiskitOptimizationError: If the problem contains non-binary variables. QiskitOptimizationError: If the problem contains constraints. """ # Ensure the Encoding object is not already used if self._frozen: raise QiskitOptimizationError( "Cannot reuse an Encoding object that has already been used. " "Please create a new Encoding object and call encode() on it." ) # Check for non-binary variables if problem.get_num_vars() > problem.get_num_binary_vars(): raise QiskitOptimizationError( "All variables must be binary. " "Please convert integer variables to binary variables using the" "``QuadraticProgramToQubo`` converter. " "Continuous variables are not supported by the QRAO algorithm." ) # Check for constraints if problem.linear_constraints or problem.quadratic_constraints: raise QiskitOptimizationError( "The problem cannot contain constraints. " "Please convert constraints to penalty terms of the objective function using the " "``QuadraticProgramToQubo`` converter." ) num_vars = problem.get_num_vars() # Generate the coefficients of the Hamiltonian offset, linear, quad = self._generate_ising_coefficients(problem) # Find the partition of the variables into groups variable_partition = self._find_variable_partition(quad) # Add variables and generate the Hamiltonian for _, v in sorted(variable_partition.items()): self._add_variables(sorted(v)) for i in range(num_vars): w = linear[i] if w != 0: self._add_term(w, i) for i in range(num_vars): for j in range(num_vars): w = quad[i, j] if w != 0: self._add_term(w, i, j) self._offset = offset self._problem = problem self.freeze()
[docs] def state_preparation_circuit(self, x: list[int]) -> QuantumCircuit: """ Generate a circuit that prepares the state corresponding to the given binary string. Args: x: A list of binary values to be encoded into the state. Returns: A QuantumCircuit that prepares the state corresponding to the given binary string. """ return _qrac_state_prep_multi_qubit(x, self.q2vars, self.max_vars_per_qubit)