r"""

The $Z_2$ symmetry converter identifies symmetries from the problem hamiltonian and uses them to

provide a tapered - more efficient - representation of operators as Paulis for this problem. For each

identified symmetry, one qubit can be eliminated in the Pauli representation at the cost of having to

test two symmetry sectors (for the two possible eigenvalues - tapering values - of the symmetry).

In certain problems such as the finding of the main operator's ground state, one can a priori

identify the symmetry sector of the solution and thus effectively reduce the computational overhead.

The following attributes can be read and updated once the ``Z2Symmetries`` object has been

constructed.

Attributes:

tapering_values (list[int] or None): Values determining the sector.

tol (float): The tolerance threshold for ignoring real and complex parts of a coefficient.

References:

[1]: Bravyi, S., et al, "Tapering off qubits to simulate fermionic Hamiltonians"

`arXiv:1701.08213 <https://arxiv.org/abs/1701.08213>`__

"""

def __init__(

self,

symmetries: Iterable[Pauli],

sq_paulis: Iterable[Pauli],

sq_list: Iterable[int],

tapering_values: Iterable[int] | None = None,

*,

tol: float = 1e-14,

):

r"""

Args:

symmetries: Object representing the list of $Z_2$ symmetries. These correspond to

the generators of the symmetry group $\langle \tau_1, \tau_2\dots \rangle>$.

sq_paulis: Object representing the list of single-qubit Pauli $\sigma^x_{q(i)}$

anti-commuting with the symmetry $\tau_i$ and commuting with all the other symmetries

$\tau_{j\neq i}$. These operators are used to construct the unitary Clifford operators.

sq_list: The list of indices $q(i)$ of the single-qubit Pauli operators used to build the

Clifford operators.

tapering_values: List of eigenvalues determining the symmetry sector for each symmetry.

tol: Tolerance threshold for ignoring real and complex parts of a coefficient.

Raises:

QiskitError: Invalid paulis. The lists of symmetries, single-qubit paulis support paulis

and tapering values must be of equal length. This length is the number of applied

symmetries and translates directly to the number of eliminated qubits.

"""

symmetries = list(symmetries)

sq_paulis = list(sq_paulis)

sq_list = list(sq_list)

tapering_values = None if tapering_values is None else list(tapering_values)

if len(symmetries) != len(sq_paulis):

raise QiskitError(

f"The number of Z2 symmetries, {len(symmetries)}, has to match the number \

of single-qubit pauli operators, {len(sq_paulis)}."

)

if len(sq_paulis) != len(sq_list):

raise QiskitError(

f"The number of single-qubit pauli operators, {len(sq_paulis)}, has to match the length \

of the of single-qubit list, {len(sq_list)}."

)

if tapering_values is not None:

if len(sq_list) != len(tapering_values):

raise QiskitError(

f"The length of the single-qubit list, {len(sq_list)}, must match the length of the \

tapering values, {len(tapering_values)} ."

)

self._symmetries = symmetries

self._sq_paulis = sq_paulis

self._sq_list = sq_list

self.tapering_values = tapering_values

self.tol = tol

@property

def symmetries(self) -> list[Pauli]:

"""Return symmetries."""

return self._symmetries

@property

def sq_paulis(self) -> list[Pauli]:

"""Return sq paulis."""

return self._sq_paulis

@property

def cliffords(self) -> list[SparsePauliOp]:

"""

Get clifford operators, built based on symmetries and single-qubit X.

Returns:

A list of unitaries used to diagonalize the Hamiltonian.

"""

cliffords = [

(SparsePauliOp(pauli_symm) + SparsePauliOp(sq_pauli)) / math.sqrt(2)

for pauli_symm, sq_pauli in zip(self._symmetries, self._sq_paulis)

]

return cliffords

@property

def sq_list(self) -> list[int]:

"""Return sq list."""

return self._sq_list

@property

def settings(self) -> dict:

"""Return operator settings."""

return {

"symmetries": self._symmetries,

"sq_paulis": self._sq_paulis,

"sq_list": self._sq_list,

"tapering_values": self.tapering_values,

}

def __str__(self):

ret = ["Z2 symmetries:"]

ret.append("Symmetries:")

for symmetry in self._symmetries:

ret.append(symmetry.to_label())

ret.append("Single-Qubit Pauli X:")

for x in self._sq_paulis:

ret.append(x.to_label())

ret.append("Cliffords:")

for c in self.cliffords:

ret.append(str(c))

ret.append("Qubit index:")

ret.append(str(self._sq_list))

ret.append("Tapering values:")

if self.tapering_values is None:

possible_values = [

str(list(coeff)) for coeff in itertools.product([1, -1], repeat=len(self._sq_list))

]

possible_values = ", ".join(x for x in possible_values)

ret.append(" - Possible values: " + possible_values)

else:

ret.append(str(self.tapering_values))

ret = "\n".join(ret)

return ret

def is_empty(self) -> bool:

"""

Check the z2_symmetries is empty or not.

Returns:

Empty or not.

"""

return len(self._symmetries) == 0 or len(self._sq_paulis) == 0 or len(self._sq_list) == 0

@classmethod

def find_z2_symmetries(cls, operator: SparsePauliOp) -> Z2Symmetries:

"""

Finds Z2 Pauli-type symmetries of a :class:`.SparsePauliOp`.

Returns:

A ``Z2Symmetries`` instance.

"""

pauli_symmetries = []

sq_paulis = []

sq_list = []

stacked_paulis = []

test_idx = {

"X_or_I": [(0, 0), (1, 0)],

"Y_or_I": [(0, 0), (1, 1)],

"Z_or_I": [(0, 0), (0, 1)],

}

test_row = {

"Z_or_I": [(1, 0), (1, 1)],

"X_or_I": [(0, 1), (1, 1)],

"Y_or_I": [(0, 1), (1, 0)],

}

pauli_bool = {

"Z_or_I": [False, True],

"X_or_I": [True, False],

"Y_or_I": [True, True],

}

if _sparse_pauli_op_is_zero(operator):

return cls([], [], [], None)

for pauli in iter(operator):

stacked_paulis.append(

np.concatenate((pauli.paulis.x[0], pauli.paulis.z[0]), axis=0).astype(int)

)

stacked_matrix = np.stack(stacked_paulis)

symmetries = _kernel_f2(stacked_matrix)

if not symmetries:

return cls([], [], [], None)

stacked_symmetries = np.stack(symmetries)

symm_shape = stacked_symmetries.shape

half_symm_shape = symm_shape[1] // 2

stacked_symm_del = [

np.delete(stacked_symmetries, row, axis=0) for row in range(symm_shape[0])

]

def _test_symmetry_row_col(row: int, col: int, idx_test: list, row_test: list) -> bool:

"""

Utility method that determines how to build the list of single-qubit Pauli X operators and

the list of corresponding qubit indices from the stacked symmetries.

This method is successively applied to Z type, X type and Y type symmetries (in this order)

to build the letter at position (col) of the Pauli word corresponding to the symmetry at

position (row).

Args:

row (int): Index of the symmetry for which the single-qubit Pauli X operator is being

built.

col (int): Index of the letter in the Pauli word corresponding to the single-qubit Pauli

X operator.

idx_test (list): List of possibilities for the stacked symmetries at all other rows

than row.

row_test (list): List of possibilities for the stacked symmetries at row.

Returns:

Whether or not this symmetry type should be used to build this letter of this

single-qubit Pauli X operator.

"""

stacked_symm_idx_tests = np.array(

[

(

stacked_symm_del[row][symm_idx, col],

stacked_symm_del[row][symm_idx, col + half_symm_shape],

)

in idx_test

for symm_idx in range(symm_shape[0] - 1)

]

)

stacked_symm_row_test = (

stacked_symmetries[row, col],

stacked_symmetries[row, col + half_symm_shape],

) in row_test

return bool(np.all(stacked_symm_idx_tests)) and stacked_symm_row_test

for row in range(symm_shape[0]):

pauli_symmetries.append(

Pauli(

(

stacked_symmetries[row, :half_symm_shape],

stacked_symmetries[row, half_symm_shape:],

)

)

)

# Try all cases for the symmetries other than row: Z or I, Y or I, X or I on col qubit.

# One test will return true.

# Build the single-qubit Pauli accordingly.

# Build the index list accordingly.

for col in range(half_symm_shape):

for key in ("Z_or_I", "X_or_I", "Y_or_I"):

current_test_result = _test_symmetry_row_col(

row, col, test_idx[key], test_row[key]

)

if current_test_result:

sq_paulis.append(

Pauli((np.zeros(half_symm_shape), np.zeros(half_symm_shape)))

)

sq_paulis[row].z[col] = pauli_bool[key][0]

sq_paulis[row].x[col] = pauli_bool[key][1]

sq_list.append(col)

break

if current_test_result:

# We break out of the loop over columns only when one valid test is identified.

break

return cls(pauli_symmetries, sq_paulis, sq_list, None)

def convert_clifford(self, operator: SparsePauliOp) -> SparsePauliOp:

"""This method operates the first part of the tapering.

It converts the operator by composing it with the clifford unitaries defined in the current

symmetry.

Args:

operator: The to-be-tapered operator.

Returns:

``SparsePauliOp`` corresponding to the converted operator.

"""

if not self.is_empty() and not _sparse_pauli_op_is_zero(operator):

# If the operator is zero then we can skip the following.

for clifford in self.cliffords:

operator = cast(SparsePauliOp, clifford @ operator @ clifford)

operator = operator.simplify(atol=0.0)

return operator

def taper_clifford(self, operator: SparsePauliOp) -> Union[SparsePauliOp, list[SparsePauliOp]]:

"""Operate the second part of the tapering.

This function assumes that the input operators have already been transformed using

:meth:`convert_clifford`. The redundant qubits due to the symmetries are dropped and

replaced by their two possible eigenvalues.

Args:

operator: Partially tapered operator resulting from a call to :meth:`convert_clifford`.

Returns:

If tapering_values is None: [:class:`SparsePauliOp`]; otherwise, :class:`SparsePauliOp`.

"""

tapered_ops: Union[SparsePauliOp, list[SparsePauliOp]]

if self.is_empty():

tapered_ops = operator

else:

# If the operator is zero we still need to taper the operator to reduce its size i.e. the

# number of qubits so for example 0*"IIII" could taper to 0*"II" when symmetries remove

# two qubits.

if self.tapering_values is None:

tapered_ops = [

self._taper(operator, list(coeff))

for coeff in itertools.product([1, -1], repeat=len(self._sq_list))

]

else:

tapered_ops = self._taper(operator, self.tapering_values)

return tapered_ops

def taper(self, operator: SparsePauliOp) -> Union[SparsePauliOp, list[SparsePauliOp]]:

"""

Taper an operator based on the z2_symmetries info and sector defined by `tapering_values`.

Returns operator if the symmetry object is empty.

The tapering is a two-step algorithm which first converts the operator into a

:class:`SparsePauliOp` with same eigenvalues but where some qubits are only acted upon

with the Pauli operators I or X.

The number M of these redundant qubits is equal to the number M of identified symmetries.

The second step of the reduction consists in replacing these qubits with the possible

eigenvalues of the corresponding Pauli X, giving 2^M new operators with M less qubits.

If an eigenvalue sector was previously identified for the solution, then this reduces to

1 new operator with M less qubits.

Args:

operator: The to-be-tapered operator.

Returns:

If tapering_values is None: [:class:`SparsePauliOp`]; otherwise, :class:`SparsePauliOp`.

"""

converted_ops = self.convert_clifford(operator)

tapered_ops = self.taper_clifford(converted_ops)

return tapered_ops

def _taper(self, op: SparsePauliOp, curr_tapering_values: list[int]) -> SparsePauliOp:

pauli_list = []

for pauli_term in iter(op):

coeff_out = pauli_term.coeffs[0]

for idx, qubit_idx in enumerate(self._sq_list):

if pauli_term.paulis.z[0, qubit_idx] or pauli_term.paulis.x[0, qubit_idx]:

coeff_out = curr_tapering_values[idx] * coeff_out

z_temp = np.delete(pauli_term.paulis.z[0].copy(), np.asarray(self._sq_list))

x_temp = np.delete(pauli_term.paulis.x[0].copy(), np.asarray(self._sq_list))

pauli_list.append((Pauli((z_temp, x_temp)).to_label(), coeff_out))

spo = SparsePauliOp.from_list(pauli_list).simplify(atol=0.0)

spo = spo.chop(self.tol)

return spo

def __eq__(self, other: Z2Symmetries) -> bool:

"""

Overload `==` operation to evaluate equality between Z2Symmetries.

Args:

other: The `Z2Symmetries` to compare to self.

Returns:

A bool equal to the equality of self and other.

"""

if not isinstance(other, Z2Symmetries):

return False

return (

self.symmetries == other.symmetries

and self.sq_paulis == other.sq_paulis

and self.sq_list == other.sq_list

and self.tapering_values == other.tapering_values

"""

Compute the kernel of a binary matrix on the binary finite field.

Args:

matrix_in (numpy.ndarray): Binary matrix.

Returns:

The list of kernel vectors.

"""

size = matrix_in.shape

kernel = []

matrix_in_id = np.vstack((matrix_in, np.identity(size[1])))

matrix_in_id_ech = (_row_echelon_f2(matrix_in_id.transpose())).transpose()

for col in range(size[1]):

if np.array_equal(

matrix_in_id_ech[0 : size[0], col], np.zeros(size[0])

) and not np.array_equal(matrix_in_id_ech[size[0] :, col], np.zeros(size[1])):

kernel.append(matrix_in_id_ech[size[0] :, col])

"""

Compute the row Echelon form of a binary matrix on the binary finite field.

Args:

matrix_in (numpy.ndarray): Binary matrix.

Returns:

Matrix_in in Echelon row form.

"""

size = matrix_in.shape

for i in range(size[0]):

pivot_index = 0

for j in range(size[1]):

if matrix_in[i, j] == 1:

pivot_index = j

break

for k in range(size[0]):

if k != i and matrix_in[k, pivot_index] == 1:

matrix_in[k, :] = np.mod(matrix_in[k, :] + matrix_in[i, :], 2)

matrix_out_temp = deepcopy(matrix_in)

indices = []

matrix_out = np.zeros(size)

for i in range(size[0] - 1):

if np.array_equal(matrix_out_temp[i, :], np.zeros(size[1])):

indices.append(i)

for row in np.sort(indices)[::-1]:

matrix_out_temp = np.delete(matrix_out_temp, (row), axis=0)

matrix_out[0 : size[0] - len(indices), :] = matrix_out_temp

matrix_out = matrix_out.astype(int)