Merge pull request #15 from QuantumApplicationLab/update_decomposition

Update decomposition

vqls_prototype/matrix_decomposition/matrix_decomposition.py

@@ -220,6 +220,16 @@ def decompose_matrix(
     ) -> Tuple[complex_array_type, List[complex_array_type], List[QuantumCircuit]]:
         raise NotImplementedError(f"can't decompose in {self.__class__.__name__!r}")
 
+    def update_matrix(self, new_matrix: npt.NDArray) -> None:
+        """Update the decomposition with a new matrix
+
+        Args:
+            new_matrix (npt.NDArray): new input matrix
+        """
+        self.sparse_matrix = spsp.issparse(new_matrix)
+        self._matrix, self.num_qubits = self._validate_matrix(new_matrix)
+        self._coefficients, self._matrices, self._circuits = self.decompose_matrix()
+
     def save(self, filename) -> None:
         """save the decomposition for future use
 
@@ -451,24 +461,55 @@ def decompose_matrix(
         possible_pauli_strings = self.get_possible_pauli_strings()
         for pauli_gates in tqdm(possible_pauli_strings):
             pauli_string = "".join(pauli_gates)
-            pauli_op = SparsePauliOp(pauli_string)  # Pauli(pauli_string)
-            # pauli_matrix = pauli_op.to_matrix()
-            # coef: complex_array_type = np.trace(pauli_matrix @ self.matrix)
-            # coef: complex_array_type = np.trace(np.dot(pauli_op, self.matrix))
-            if self.sparse_matrix:
-                coef: complex_array_type = (
-                    pauli_op.to_matrix(sparse=True) @ self.matrix
-                ).trace()
-            else:
-                coef: complex_array_type = np.einsum("ij,ji", pauli_op, self.matrix)  # type: ignore
-
+            coef: complex_array_type = self._get_pauli_coefficient(
+                self.matrix, pauli_string, self.sparse_matrix
+            )
             if coef * np.conj(coef) != 0:
                 self.strings.append(pauli_string)
                 coeffs.append(prefactor * coef)
                 circuits.append(self._create_circuit(pauli_string))
 
         return np.array(coeffs, dtype=np.cdouble), unit_mats, circuits
 
+    def update_matrix(self, new_matrix: npt.NDArray) -> None:
+        """Update the coefficients using a new matrix
+
+        Args:
+            new_matrix (npt.NDArray): new input matrix
+
+        """
+        prefactor = 1.0 / (2**self.num_qubits)
+        coeffs = []
+        for pauli_string in tqdm(self.strings):
+            coef = self._get_pauli_coefficient(
+                new_matrix, pauli_string, self.sparse_matrix
+            )
+            coeffs.append(prefactor * coef)
+
+        self._matrix = new_matrix
+        self.coefficients = np.array(coeffs, dtype=np.cdouble)
+
+    @staticmethod
+    def _get_pauli_coefficient(
+        matrix: npt.ArrayLike, pauli_string: str, sparse_matrix: bool
+    ) -> complex_array_type:
+        """Compute the pauli coefficient of a given pauli string
+
+        Args:
+            matrix (npt.ArrayLike): input matrix
+            pauli_string (str): a given pauli string
+            sparse_matrix (bool): if the matrix is sparse or not
+
+        Returns:
+            complex_array_type: pauli coefficient
+        """
+        pauli_op = SparsePauliOp(pauli_string)  # Pauli(pauli_string)
+        if sparse_matrix:
+            coef = (pauli_op.to_matrix(sparse=True) @ matrix).trace()
+        else:
+            coef = np.einsum("ij,ji", pauli_op, matrix)  # type: ignore
+        return coef
+
     def recompose(self) -> complex_array_type:
         """Recompose the full matrix
 

vqls_prototype/solver/vqls.py

@@ -179,6 +179,17 @@ def __init__(
         }
         self.options = self._validate_solve_options(options)
 
+        self.supported_decomposition = {
+            "pauli": PauliDecomposition,
+            "contracted_pauli": ContractedPauliDecomposition,
+            "optimized_pauli": OptimizedPauliDecomposition,
+            "symmetric": SymmetricDecomposition,
+        }
+
+        self.supported_decomposition_list = tuple(
+            v for _, v in self.supported_decomposition.items()
+        )
+
     def construct_circuit(  # pylint: disable=too-many-branches
         self,
         matrix: Union[np.ndarray, QuantumCircuit, List],
@@ -224,13 +235,16 @@ def construct_circuit(  # pylint: disable=too-many-branches
         else:
             raise ValueError("Format of the input vector not recognized")
 
-        # general numpy matrix
+        # Reuse the matrix if we reinit with a different rhs
         if (self.options["reuse_matrix"] is True) and (
             self.matrix_circuits is not None
         ):
             print("\t VQLS : Reusing matrix decomposition for new RHS")
 
+        # recreate a decomposition for the matrix
         else:
+
+            # general np array
             if isinstance(matrix, np.ndarray):
                 # ensure the matrix is double
                 matrix = matrix.astype("float64")
@@ -243,14 +257,15 @@ def construct_circuit(  # pylint: disable=too-many-branches
                         + ". Matrix dimension: "
                         + str(matrix.shape[0])
                     )
-                decomposition = {
-                    "pauli": PauliDecomposition,
-                    "contracted_pauli": ContractedPauliDecomposition,
-                    "optimized_pauli": OptimizedPauliDecomposition,
-                    "symmetric": SymmetricDecomposition,
-                }[self.options["matrix_decomposition"]]
+                decomposition = self.supported_decomposition[
+                    self.options["matrix_decomposition"]
+                ]
                 self.matrix_circuits = decomposition(matrix=matrix)
 
+            # a pregenerated decomposition
+            elif isinstance(matrix, self.supported_decomposition_list):
+                self.matrix_circuits = matrix  # type: ignore[assignment]
+
             # a single circuit
             elif isinstance(matrix, QuantumCircuit):
                 if matrix.num_qubits != self.vector_circuit.num_qubits:
@@ -268,7 +283,9 @@ def construct_circuit(  # pylint: disable=too-many-branches
                 )
 
             else:
-                raise ValueError("Format of the input matrix not recognized")
+                raise ValueError(
+                    "Format of the input matrix not recognized", type(matrix)
+                )
 
         # create only the circuit for <psi|psi> =  <0|V A_n ^* A_m V|0>
         # with n != m as the diagonal terms (n==m) always give a proba of 1.0