updated from_matrix method of PauliWordOp. Fixed bugs, faster code an…

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> <j|
+
+    Args:
+        i (int): ket of projector
+        j (int): bra of projector
+        n_qubits (int): number of qubits
+        binary_vec (optional): bool array of all bitstrings on n_qubits
+        return_operator (bool): whether to return PauliWordOp or a the symplectic matrix and coeff_vec
+
     """
     if n_qubits > 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:

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