diff --git a/symmer/operators/base.py b/symmer/operators/base.py index e74a12fa..0983e207 100644 --- a/symmer/operators/base.py +++ b/symmer/operators/base.py @@ -9,7 +9,7 @@ from numbers import Number import multiprocessing as mp from cached_property import cached_property -from scipy.sparse import csr_matrix, csc_matrix, coo_matrix +from scipy.sparse import csr_matrix, csc_matrix, coo_matrix, dok_matrix from symmer.operators.utils import * from openfermion import QubitOperator, count_qubits import matplotlib.pyplot as plt @@ -171,23 +171,33 @@ def _from_matrix_full_basis(cls, XZ_block = (((int_list[:, None] & (1 << np.arange(2 * n_qubits))[::-1])) > 0).astype(int) op_basis = cls(XZ_block, np.ones(XZ_block.shape[0])) else: - op_basis = operator_basis + op_basis = operator_basis.cleanup().copy() + op_basis.coeff_vec = np.ones(op_basis.coeff_vec.shape) denominator = 2 ** n_qubits decomposition = cls.empty(n_qubits) for op in tqdm(op_basis, desc='Building operator via full basis', total=op_basis.n_terms): if isinstance(matrix, np.ndarray): - const = np.einsum( - 'ij,ij->', - op.to_sparse_matrix.toarray(), - matrix, - optimize=True - ) / denominator + ### dense operation! + # const = np.einsum( + # 'ij,ij->', + # op.to_sparse_matrix.toarray(), + # matrix, + # optimize=True + # ) / denominator + const = np.multiply(op.to_sparse_matrix.toarray(), matrix).sum() / denominator else: + ### sparse operation! const = (op.to_sparse_matrix.multiply(matrix)).sum() / denominator + decomposition += op.multiply_by_constant(const) operator_out = decomposition.cleanup() + + # fix ZX Y phases generated! + Y_sign = (operator_out.Y_count % 2 * -2) + 1 + operator_out.coeff_vec = operator_out.coeff_vec * Y_sign + if operator_basis is not None: if not np.all(operator_out.to_sparse_matrix.toarray() == matrix): warnings.warn('Basis not sufficiently expressive, output operator projected onto basis supplied.') @@ -201,25 +211,49 @@ def _from_matrix_projector(cls, ) -> "PauliwordOp": """ """ + assert n_qubits <= 32, 'cannot decompose matrices above 32 qubits' + if isinstance(matrix, np.ndarray): row, col = np.where(matrix) elif isinstance(matrix, (csr_matrix, csc_matrix, coo_matrix)): row, col = matrix.nonzero() else: raise ValueError('Unrecognised matrix type, must be one of np.array or sp.sparse.csr_matrix') - - binary_vec = ( - ( - np.arange(2 ** n_qubits).reshape([-1, 1]) & - (1 << np.arange(n_qubits))[::-1] - ) > 0 - ).astype(bool) - P_out = cls.empty(n_qubits) - for i,j in tqdm(zip(row, col), desc='Building operator via projectors', total=len(row)): - ij_op = get_ij_operator(i,j,n_qubits,binary_vec=binary_vec) - P_out += ij_op * matrix[i,j] + sym_operator = dok_matrix((4 ** n_qubits, 2 * n_qubits), + dtype=bool) + + coeff_operator = dok_matrix((4 ** n_qubits, 1), + dtype=complex) + binary_vec = ( + ( + np.arange(2 ** n_qubits).reshape([-1, 1]) & + (1 << np.arange(n_qubits))[::-1] + ) > 0).astype(bool) + + binary_convert = 1 << np.arange(2 * n_qubits)[::-1] + # P_out = cls.empty(n_qubits) + for i, j in tqdm(zip(row, col), desc='Building operator via projectors', total=len(row)): + ij_symp_matrix, proj_coeffs = get_ij_operator(i, j, + n_qubits, + binary_vec=binary_vec, + return_operator=False) + + ### find location in symp matrix + int_list = ij_symp_matrix @ binary_convert # (1 << np.arange(ij_symp_matrix.shape[1])[::-1]) + + # populate sparse mats + sym_operator[int_list, :] = ij_symp_matrix + coeff_operator[int_list] += proj_coeffs.reshape(-1, 1) * matrix[i, j] + + ### only keep nonzero coeffs! (skips expensive cleanup) + nonzero = coeff_operator.nonzero()[0] + P_out = PauliwordOp(sym_operator[nonzero, :].toarray(), + coeff_operator[nonzero].toarray()[:, 0]) + + # P_out = PauliwordOp(sym_operator.toarray(), + # coeff_operator.toarray()[:,0]).cleanup() return P_out @classmethod @@ -1617,41 +1651,60 @@ def get_PauliwordOp_projector(projector: Union[str, List[str], np.array]) -> "Pa return projector -def get_ij_operator(i,j,n_qubits,binary_vec=None): +def get_ij_operator(i:int, j:int, n_qubits:int, + binary_vec:np.ndarray=None, + return_operator:bool=True) -> Union["PauliwordOp", Tuple[np.ndarray, np.ndarray]]: """ + Get the Pauli operator for the projector: |i> 30: raise ValueError('Too many qubits, might run into memory limitations.') if binary_vec is None: binary_vec = ( - ((np.arange(2 ** n_qubits).reshape([-1, 1]) & - (1 << np.arange(n_qubits))[::-1])) > 0 + ((np.arange(2 ** n_qubits).reshape([-1, 1]) & + (1 << np.arange(n_qubits))[::-1])) > 0 ).astype(bool) - left = np.array([int(i) for i in np.binary_repr(i, width=n_qubits)]).astype(bool) + left = np.array([int(i) for i in np.binary_repr(i, width=n_qubits)]).astype(bool) right = np.array([int(i) for i in np.binary_repr(j, width=n_qubits)]).astype(bool) - AND = left & right # AND where -1 sign - XZX_sign_flips = (-1) ** np.sum(AND & binary_vec, axis=1) # XZX = -X multiplications - + AND = left & right # AND where -1 sign + XZX_sign_flips = (-1) ** np.sum(AND & binary_vec, axis=1) # XZX = -X multiplications + if i != j: - XOR = left ^ right # XOR where +-i phase + XOR = left ^ right # XOR where +-i phase XZ_mult = left & binary_vec ZX_mult = binary_vec & right - XZ_phase = (-1j) ** np.sum(XZ_mult & ~ZX_mult, axis=1) # XZ=-iY multiplications - ZX_phase = (+1j) ** np.sum(ZX_mult & ~XZ_mult, axis=1) # ZX=+iY multiplications + XZ_phase = (-1j) ** np.sum(XZ_mult & ~ZX_mult, axis=1) # XZ=-iY multiplications + ZX_phase = (+1j) ** np.sum(ZX_mult & ~XZ_mult, axis=1) # ZX=+iY multiplications phase_mod = XZX_sign_flips * XZ_phase * ZX_phase - - ij_symp_matrix = np.hstack([np.tile(XOR, [2**n_qubits, 1]), binary_vec]) - ij_operator= PauliwordOp(ij_symp_matrix, phase_mod/2**n_qubits) + + ij_symp_matrix = np.hstack([np.tile(XOR, [2 ** n_qubits, 1]), binary_vec]) + coeffs = phase_mod / 2 ** n_qubits + + if return_operator: + ij_operator = PauliwordOp(ij_symp_matrix, phase_mod / 2 ** n_qubits) + return ij_operator else: ij_symp_matrix = np.hstack([np.zeros_like(binary_vec), binary_vec]) - ij_operator= PauliwordOp(ij_symp_matrix, XZX_sign_flips/2**n_qubits) - - return ij_operator + coeffs = XZX_sign_flips / 2 ** n_qubits + + if return_operator: + ij_operator = PauliwordOp(ij_symp_matrix, XZX_sign_flips / 2 ** n_qubits) + return ij_operator + + return ij_symp_matrix, coeffs def single_term_expval(P_op: PauliwordOp, psi: QuantumState) -> float: diff --git a/tests/test_operators/test_base.py b/tests/test_operators/test_base.py index e76a49ef..f3d0369d 100644 --- a/tests/test_operators/test_base.py +++ b/tests/test_operators/test_base.py @@ -2,6 +2,7 @@ import numpy as np from symmer.operators import PauliwordOp from functools import reduce +from scipy.sparse import rand, csr_matrix P_matrices ={ 'X': np.array([[0, 1], @@ -260,25 +261,131 @@ def test_to_dictionary( ).to_dictionary == pauli_dict -def test_from_matrix( - pauli_list_1, - symp_matrix_1, - coeff_vec_1 - ): - n_qubits = len(pauli_list_1[0]) - - # build matrix - mat = np.zeros((2**n_qubits, 2**n_qubits), dtype=complex) - for index, p_op in enumerate(pauli_list_1): - coeff = coeff_vec_1[index] - p_op_mat = reduce(np.kron, [P_matrices[sig] for sig in p_op]) - mat+= coeff*p_op_mat - # generate Pop from matrix - PauliOp_from_matrix = PauliwordOp.from_matrix(mat) +def test_from_matrix_projector_dense(): + + for n_qubits in range(1,5): + mat = np.random.random((2**n_qubits,2**n_qubits)) + 1j*np.random.random((2**n_qubits,2**n_qubits)) + PauliOp_from_matrix = PauliwordOp.from_matrix(mat, + strategy='projector') + assert np.allclose(PauliOp_from_matrix.to_sparse_matrix.toarray(), + mat) + + +def test_from_matrix_full_basis_dense(): + + for n_qubits in range(1,5): + mat = np.random.random((2**n_qubits,2**n_qubits)) + 1j*np.random.random((2**n_qubits,2**n_qubits)) + PauliOp_from_matrix = PauliwordOp.from_matrix(mat, + strategy='full_basis') + + assert np.allclose(PauliOp_from_matrix.to_sparse_matrix.toarray(), + mat) + + +def test_from_matrix_defined_basis_dense(): + + + for n_qubits in range(2,6): + n_terms = (4**n_qubits)//2 + op_basis = PauliwordOp.random(n_qubits, n_terms) + PauliOp_from_matrix = PauliwordOp.from_matrix(op_basis.to_sparse_matrix.toarray(), + strategy='full_basis', + operator_basis=op_basis) + assert np.allclose(op_basis.to_sparse_matrix.toarray(), + PauliOp_from_matrix.to_sparse_matrix.toarray()) + + + op_basis = PauliwordOp.from_dictionary({'XX':1, + 'ZZ':2, + 'YY':2, + 'YI':-1, + 'YY':2, + 'ZX':2}) + + mat = np.array([[ 2.+0.j, 0.+0.j, 0.+1.j, -1.+0.j], + [ 0.+0.j, -2.+0.j, 3.+0.j, 0.+1.j], + [ 0.-1.j, 3.+0.j, -2.+0.j, 0.+0.j], + [-1.+0.j, 0.-1.j, 0.+0.j, 2.+0.j]]) + PauliOp_from_matrix = PauliwordOp.from_matrix(mat, strategy='full_basis', operator_basis=op_basis) + assert np.allclose(PauliOp_from_matrix.to_sparse_matrix.toarray(), mat) + + +def test_from_matrix_projector_sparse(): + density = 0.8 + for n_qubits in range(1,5): + dim = 2 ** n_qubits + mat = rand(dim, dim, + density=density, + format='csr', + dtype=complex) + PauliOp_from_matrix = PauliwordOp.from_matrix(mat, + strategy='projector') + assert np.allclose(PauliOp_from_matrix.to_sparse_matrix.toarray(), + mat.toarray()) + + +def test_from_matrix_full_basis_sparse(): + + density = 0.8 + for n_qubits in range(1,5): + dim = 2 ** n_qubits + mat = rand(dim, dim, + density=density, + format='csr', + dtype=complex) + PauliOp_from_matrix = PauliwordOp.from_matrix(mat, + strategy='full_basis') + + assert np.allclose(PauliOp_from_matrix.to_sparse_matrix.toarray(), + mat.toarray()) + + +def test_from_matrix_defined_basis_sparse(): + + op_basis = PauliwordOp.from_dictionary({ + 'II': 1, + 'IZ': 1, + 'ZI': 1, + 'ZZ': 1, + 'IX': 1, + 'IY': 1, + 'ZX': 1, + 'ZY': 1, + 'XI': 1, + 'XZ': 1, + 'YI': 1, + 'YZ': 1, + 'XX': 1, + 'XY': 1, + 'YX': 1, + 'YY': 1 }) + + mat = np.array([[1, 0, 0, 0], + [1, 0, 0, 0], + [1, 0, -1j, 0], + [1, 0, 0, 0]]) + sparse_mat = csr_matrix(mat) + PauliOp_from_matrix = PauliwordOp.from_matrix(sparse_mat, + strategy='full_basis', + operator_basis=op_basis) + assert np.allclose(PauliOp_from_matrix.to_sparse_matrix.toarray(), mat) + + +def test_from_matrix_incomplete_op_basis(): + """ + Test to see if warning thrown if supplied basis is not enough. + Returns: + + """ + op_basis = PauliwordOp.from_dictionary({'XX': 1}) + mat = np.array([[2. + 0.j, 0. + 0.j, 0. + 1.j, -1. + 0.j], + [0. + 0.j, -2. + 0.j, 3. + 0.j, 0. + 1.j], + [0. - 1.j, 3. + 0.j, -2. + 0.j, 0. + 0.j], + [-1. + 0.j, 0. - 1.j, 0. + 0.j, 2. + 0.j]]) + + with pytest.warns(UserWarning): + PauliOp_from_matrix = PauliwordOp.from_matrix(mat, operator_basis=op_basis) - pauli_dict = dict(zip(pauli_list_1, coeff_vec_1)) - PauliOp_from_dict = PauliwordOp.from_dictionary(pauli_dict) - assert PauliOp_from_matrix == PauliOp_from_dict def test_from_matrix_to_matrix(): n_qubits = 3