diff --git a/h4_experiment.py b/h4_experiment.py index 08cb1d6..07dac9b 100644 --- a/h4_experiment.py +++ b/h4_experiment.py @@ -1,7 +1,9 @@ +from typing import Dict, List, Tuple import os import pickle from tqdm import tqdm import numpy as np +import scipy from multiprocessing import Pool import datetime @@ -17,8 +19,15 @@ from sampler import local_dists_optimal -def _convert_spin_sector(ham_f): - "convert |↑↓↑↓...> representation to |↑↑...↓↓...> representation" +def _convert_spin_sector(ham_f: FermionOperator) -> FermionOperator: + """convert |↑↓↑↓...> representation to |↑↑...↓↓...> representation. + + Args: + ham_f (FermionOperator): Second-quantized molecular Hamiltonian in |↑↓↑↓...> representation. + + Returns: + FermionOperator: Second-quantized molecular Hamiltonian in |↑↑...↓↓...> representation. + """ n = count_qubits(ham_f) n_half = int(n / 2) index_map = {i: int(i / 2) if i % 2 == 0 else n_half + int((i - 1) / 2) for i in range(n)} @@ -27,7 +36,15 @@ def _convert_spin_sector(ham_f): return ham_convert -def _get_h4_hamiltonian(length): +def _get_h4_hamiltonian(length: float) -> QubitOperator: + """Generate H4 Hamiltonian from given interatomic distance. + + Args: + length (float): Interatomic distance (Å). + + Returns: + QubitOperator: H4 Hamiltonian. + """ geom = [ ("H", (-length / 2 * 3, 0, 0)), ("H", (-length / 2, 0, 0)), @@ -52,7 +69,13 @@ def _get_h4_hamiltonian(length): return ham -def simulate_energy_vs_interatomic_distance(): +def simulate_energy_vs_interatomic_distance() -> Tuple[float, float, float]: + """Calculate energies with various interatomic distances. + + Returns: + Tuple[float, float, float]: Rigorous energies, and energies estimated through the ground state and CISD state. + """ + energies_rigorous = [] energies_2n1p_gs = [] energies_2n1p_cisd = [] @@ -79,7 +102,17 @@ def simulate_energy_vs_interatomic_distance(): return energies_rigorous, energies_2n1p_gs, energies_2n1p_cisd -def _estimate_qse_matrix_elements_noise_simulator(beta_eff, params): +def _estimate_qse_matrix_elements_noise_simulator(beta_eff: np.ndarray, params: Dict) -> Tuple[np.ndarray, np.ndarray]: + """Estimate operator matrices of QSE through LBCS using noise simulator. + + Args: + beta_eff (np.ndarray): Weighted bias of measurement bases within LBCS. + params (Dict): Parameters for execution. + + Returns: + Tuple[np.ndarray, np.ndarray]: Estimated matrix elements \tilde{H}_ij = <\psi| O_i^\dag H O_j |\psi> and \tilde{S}_ij = <\psi| O_i^\dag O_j |\psi> + """ + def estimate_qse_matrix_element_noise_simulator(i, j, op_mat, result_mat): if j >= i: for term, coef in op_mat[i, j].terms.items(): @@ -116,7 +149,15 @@ def get_n_shots_per_pauli(term, beta, shots_total): return h_eff, s_mtrc -def _iterative_run_lbcs(params): +def _iterative_run_lbcs(params: Dict) -> Dict: + """Apply adaptive measurement optimization (LBCS) + + Args: + params (Dict): Parameters for execution. + + Returns: + Dict: QSE results in each iteration + """ subspace_expansion = params["subspace_expansion"] hamiltonian = params["molecule"].hamiltonian n_qubit = subspace_expansion.n_qubit @@ -171,7 +212,20 @@ def _iterative_run_lbcs(params): return dict_iter_lbcs -def _get_pauli_mat_dict(file_name, n_qubits, h_ij, load=True): +def _get_pauli_mat_dict( + file_name: str, n_qubits: int, h_ij: Dict[Tuple, QubitOperator], load=True +) -> Dict[Tuple, scipy.sparse.csc.csc_matrix]: + """Generate dictionary of sparse matrix. + + Args: + file_name (str): Filename of sparse matrix dictionary to loda/save. + n_qubits (int): Number of qubits. + h_ij (Dict[Tuple, QubitOperator]): Operator matrix + load (bool, optional): If True, precomputed sparse matrix dictionary is used. Defaults to True. + + Returns: + Dict[Tuple, scipy.sparse]: _description_ + """ if load and file_name in os.listdir("h4_experiment_results"): with open(f"h4_experiment_results/{file_name}", "rb") as f: pauli_mat_dict = pickle.load(f) @@ -185,7 +239,12 @@ def _get_pauli_mat_dict(file_name, n_qubits, h_ij, load=True): return pauli_mat_dict -def simulate_qse_convergence(): +def simulate_qse_convergence() -> Tuple[float, float, List[float]]: + """Benchmark convergence of adaptive update of measurement strategy + + Returns: + Tuple[float, float, List[float]]: rigorous and CISD energies, and energies for iterations. + """ params_qse = { "molecule": "H4", "n_qubits": 8, diff --git a/reproduce_experiment.ipynb b/reproduce_experiment.ipynb index 439362f..edbd918 100644 --- a/reproduce_experiment.ipynb +++ b/reproduce_experiment.ipynb @@ -30,6 +30,22 @@ "# Table 1, 3 and Fig. 2" ] }, + { + "cell_type": "markdown", + "id": "e651db91", + "metadata": {}, + "source": [ + "Description of Parameters:\n", + "* n_trial [int]: Number of trials to gather statistics for mean absolute error and standard deviation.\n", + "* n_lev[\"auto\" or int]: Regularization parameter for the general eigenvalue problem. The maximum value is the dimension of the subspace. If \"auto\" is selected, n_lev is chosen to minimize the absolute error of estimation.\n", + "* subspace [\"1n\" or \"2n1p\"]: Choosing \"1n\" yields a subspace $S = Span\\{a_i |\\psi\\rangle \\}$, while \"2n1p\" yields a subspace $S = Span\\{a_i a_j^\\dagger a_k |\\psi\\rangle \\}$.\n", + "* spin_supspace [\"up\", \"down\", or \"all\"]: Selecting \"up(down)\" yields a supspace $S = Span\\{O_{i,\\sigma=\\uparrow(\\downarrow) } |\\psi\\rangle \\} $, while \"all\" yields a subspace $S = Span\\{O_{i,\\sigma\\in\\{\\uparrow, \\downarrow \\}} |\\psi\\rangle \\} $.\n", + "* cpu_assigned [int]: Number of CPUs assigned for the calculation.\n", + "* verbose [0, 1, 2]: Verbosity level for the calculation progress. To minimize output, set it to 0.\n", + "* load [bool]: If True, precomputed values are utilized for some calculations.\n", + "* write_result_matrix [bool]: If True, $\\tilde{H}_{ij}$ and $\\tilde{S}_{ij}$ are saved after the calculation.\"" + ] + }, { "cell_type": "code", "execution_count": 17, @@ -40,8 +56,8 @@ "outputs": [], "source": [ "params_fixed = {\n", - " \"n_trial\": 10,\n", - " \"n_lev\": \"auto\",\n", + " \"n_trial\": 10, \n", + " \"n_lev\": \"auto\", # parameter \n", " \"subspace\": \"1n\",\n", " \"spin_supspace\": \"up\",\n", " \"cpu_assigned\": 10,\n", diff --git a/sampler.py b/sampler.py index efa6e4a..ffd36c7 100644 --- a/sampler.py +++ b/sampler.py @@ -1,4 +1,4 @@ -from typing import Iterable, List, Optional +from typing import Dict, Iterable, List, Optional import numpy as np import pandas as pd @@ -35,19 +35,22 @@ def generate_random_measurement_axis( self, lbcs_beta: Optional[Iterable[Iterable]] = None, ogm_meas_set: Optional[Iterable[Iterable]] = None, - ): - """ - Creates a list of random measurement axis. + ) -> np.ndarray: + """ Creates a list of random measurement axis. For 3-qubit system, this looks like, e.g., [1, 3, 2], which tells you to measure in [X, Z, Y] basis. - return: - List[n_qubit] + Args: + lbcs_beta (Optional[Iterable[Iterable]], optional): Weighted bias of measurement bases for LBCS. Defaults to None. + ogm_meas_set (Optional[Iterable[Iterable]], optional): Probability distribution of basis for OGM. Defaults to None. + + Returns: + np.ndarray: Optimized measurement bases. """ - assert not ((lbcs_beta is not None) and (ogm_meas_set is not None)), "you must choose either of lbcs or ogm" + assert not ((lbcs_beta is not None) and (ogm_meas_set is not None)), "you must choose either of LBCS or OGM" if lbcs_beta is not None: meas_axes = np.vstack( [np.random.multinomial(1, beta_i, size=self.Ntot).argmax(axis=1) + 1 for beta_i in lbcs_beta] @@ -60,12 +63,15 @@ def generate_random_measurement_axis( meas_axes = np.random.randint(1, 4, size=(self.Ntot, self.n_qubit)) return meas_axes - def _sample_digits(self, meas_axis, nshot_per_axis=1): - """ - returns the measurement result at meas_axis. - attributes: - meas_axis: List of axes - nshot_per_axis: number of measurement per axis + def _sample_digits(self, meas_axis: np.ndarray, nshot_per_axis=1) -> List[List[int]]: + """ Returns the measurement result at meas_axis. + + Args: + meas_axis (np.ndarray): Optimized measurement basis axis. + nshot_per_axis (int, optional): number of measurement per axis. Defaults to 1. + + Returns: + List[List[int]]: List of measurement result in same format of meas_axis. """ meas_state = QuantumState(self.n_qubit) meas_state.load(self.state) @@ -89,9 +95,26 @@ def _sample_digits(self, meas_axis, nshot_per_axis=1): return digits -def local_dists_optimal(ham, num_qubits, objective, method, β_initial=None, bitstring_HF=None): +def local_dists_optimal( + ham: QubitOperator, + num_qubits: int, + objective: str, + method: str, + β_initial: Dict = None, + bitstring_HF: str = None, +) -> np.ndarray: """Find optimal probabilities beta_{i,P} and return as dictionary - attn: qiskit ordering""" + + Args: + ham (QubitOperator): Target Hamiltonian + num_qubits (int): Number of qubits + objective (str): Objective fuction + method (str): Optimization method + bitstring_HF (str, optional): HF representation. Defaults to None. + + Returns: + np.ndarray: _description_ + """ assert objective in ["diagonal", "mixed"] assert method in ["scipy", "lagrange"] @@ -138,12 +161,13 @@ def estimate_exp( meas_axes: Iterable[Iterable] = None, samples: np.ndarray = None, ) -> float: - """ - Estimate expectation value of Observable for Basic Classical Shadow + """Estimate expectation value of Observable for Basic Classical Shadow Args: operator (QubitOperator): Observable such as Hamiltonian sampler (LocalPauliShadowSampler_core): Sampler Class + meas_axes (Iterable[Iterable], optional): Precomputed measurement axes. Defaults to None. + samples (np.ndarray, optional): Precomputed sampling results for given measurement axes. Defaults to None. Returns: float: Expectation value @@ -179,13 +203,15 @@ def estimate_exp_lbcs( meas_axes: Iterable[Iterable] = None, samples: np.ndarray = None, ) -> float: - """ - Estimate expectation value of Observable for Locally Biased Classical Shadow + """Estimate expectation value of Observable for Locally Biased Classical Shadow Args: operator (QubitOperator): Observable such as Hamiltonian sampler (LocalPauliShadowSampler_core): Sampler Class beta (Iterable[Iterable]): weighted bias + meas_axes (Iterable[Iterable], optional): Precomputed measurement axes. Defaults to None. + samples (np.ndarray, optional): Precomputed sampling results for given measurement axes. Defaults to None. + Returns: float: Expectation value @@ -231,6 +257,9 @@ def estimate_exp_ogm( sampler (LocalPauliShadowSampler_core): Sampler Class meas_dist (Iterable[Iterable]): measurement set and its distribution input is format of QubitOperator.terms + meas_axes (Iterable[Iterable], optional): Precomputed measurement axes. Defaults to None. + samples (np.ndarray, optional): Precomputed sampling results for given measurement axes. Defaults to None. + Returns: float: Expectation value @@ -287,6 +316,8 @@ def estimate_exp_derand( Args: operator (QubitOperator): Observable such as Hamiltonian sampler (LocalPauliShadowSampler_core): Sampler Class + meas_axes (Iterable[Iterable], optional): Precomputed measurement axes. Defaults to None. + samples (np.ndarray, optional): Precomputed sampling results for given measurement axes. Defaults to None. Returns: float: Expectation value diff --git a/subspace_expansion.py b/subspace_expansion.py index e595330..2adce1a 100644 --- a/subspace_expansion.py +++ b/subspace_expansion.py @@ -1,4 +1,4 @@ -from typing import Dict +from typing import Dict, Tuple, Union from multiprocessing import Pool import os import pickle @@ -55,8 +55,20 @@ def _get_qse_ops(self): terms = [FermionOperator(f"{i} {j}^ {k}") for i, j, k in itertools.product(iterator, repeat=3)] return [jordan_wigner(op) for op in terms] - def solve_qde_classically(self, molecule, state): - # preprocess (classical) subspace-expansion + def solve_qde_classically(self, molecule: MoleculeInfo, state: np.ndarray) -> Tuple[np.ndarray, float]: + """ + Performing classical preprocessing to approximate the coefficients of variational parameters. + + Args: + molecule (MoleculeInfo): Target molecule information. + state (np.ndarray): Classically simulatable state. + + Raises: + ValueError: Raised when the number of shots is too low and n_lev is too large. + + Returns: + Tuple[np.ndarray, float]: Variational coefficients and estimated energy in the classical regime. + """ terms_eff_dict = {} for i, term in enumerate(self.qse_ops): terms_eff_dict[str(term)] = get_sparse_operator(term, n_qubits=self.n_qubit) @ state @@ -97,7 +109,17 @@ def solve_qde_classically(self, molecule, state): return alpha, energy_excited - def _get_h_s_d(self, alpha_se, ham): + def _get_h_s_d(self, alpha_se: np.ndarray, ham: QubitOperator) -> Tuple[QubitOperator, QubitOperator]: + """ + Calculate dressed operators from the Hamiltonian and variational coefficients. + + Args: + alpha_se (np.ndarray): Variational coefficients for subspace expansion. + ham (QubitOperator): Molecular Hamiltonian. + + Returns: + Tuple[QubitOperator, QubitOperator]: Dressed operators H_d and S_d in Jordan-Wigner representation. + """ def qubit_op_from_term(terms): op = QubitOperator() op.terms = terms @@ -111,7 +133,17 @@ def qubit_op_from_term(terms): return H_d_jw, S_d_jw - def _get_h_s_ij(self, ham): + def _get_h_s_ij( + self, ham: QubitOperator + ) -> Tuple[Dict[Tuple[int, int], QubitOperator], Dict[Tuple[int, int], QubitOperator]]: + """Calculate operator matrices H_ij = O_i^\dag H O_j and S_ij = O_i^\dag O_j + + Args: + ham (QubitOperator): Molecular Hamiltonian. + + Returns: + Tuple[Dict[Tuple[int, int], QubitOperator], Dict[Tuple[int, int], QubitOperator]]: Operator matrices H_ij and S_ij + """ h_ij = {} s_ij = {} for i, term_i in enumerate(self.qse_ops): @@ -124,7 +156,15 @@ def _get_h_s_ij(self, ham): s_ij[i, j] = s_ij[j, i] return h_ij, s_ij - def estimate_qse_matrix_elements(self, meas_axes): + def estimate_qse_matrix_elements(self, meas_axes: np.ndarray) -> Tuple[np.ndarray, np.ndarray]: + """Estimate each element of operator matrices using measurement basis optimization subroutines + + Args: + meas_axes (np.ndarray): Optimized measurement bases + + Returns: + Tuple[np.ndarray, np.ndarray]: Estimated matrix elements \tilde{H}_ij = <\psi| O_i^\dag H O_j |\psi> and \tilde{S}_ij = <\psi| O_i^\dag O_j |\psi> + """ def estimate_qse_matrix_element(i, j, op_mat, result_mat): if j >= i: if self.params["method"] == "naive_LBCS": @@ -167,11 +207,28 @@ def estimate_qse_matrix_element(i, j, op_mat, result_mat): return H_eff, S_mtrc def _solve_nlev_regularized_gen_eig( - self, h, s, n_lev=None, threshold=1e-15, vec_norm_thresh=np.infty, return_vec=False - ): - s_vals, s_vecs = scipy.linalg.eigh( - s, - ) + self, + h: np.ndarray, + s: np.ndarray, + n_lev: int = None, + threshold=1e-15, + vec_norm_thresh=np.infty, + return_vec=False, + ) -> Union[float, Tuple[float, np.ndarray]]: + """Solve the general eigenvalue problem using regularization + + Args: + h (np.ndarray): Estimated matrix elements \tilde{H}_ij + s (np.ndarray): Estimated matrix elements \tilde{S}_ij + n_lev (int, optional): Regularization parameter (truncation size). Defaults to None. + threshold (float, optional): Regularization parameter (truncation lower bound). Defaults to 1e-15. + vec_norm_thresh (float, optional): Regularization parameter (truncation upper bound).. Defaults to np.infty. + return_vec (bool, optional): If True, return value includes the corresponding eigenvector. Defaults to False. + + Returns: + Union[float, (Tuple[float, np.ndarray]]: The lowest eigenvalue (and its eigenvector if return_vec=True) + """ + s_vals, s_vecs = scipy.linalg.eigh(s) s_vecs = s_vecs.T if n_lev is not None: good_vecs = np.array( @@ -188,14 +245,34 @@ def _solve_nlev_regularized_gen_eig( return sol[0][lowest_stable_index], sol[1][:, lowest_stable_index] @ good_vecs return sol[0][lowest_stable_index] - def _solve_regularized_gen_eig_with_best_n_lev(self, h_eff, s_mtrc, true_energy_excited, threshold, return_vec): + def _solve_regularized_gen_eig_with_best_n_lev( + self, + h: np.ndarray, + s: np.ndarray, + true_energy_excited: float, + threshold: float, + return_vec: bool, + ) -> Union[float, Tuple[float, np.ndarray]]: + """Solve the general eigenvalue problem using regularization with choosing n_lev to minimize the absolute error given a rigorous solution. + Note that this method is for benchmarking purposes and is not available when the rigorous solution is unknown. + + Args: + h (np.ndarray): Estimated matrix elements \tilde{H}_ij + s (np.ndarray): Estimated matrix elements \tilde{S}_ij + true_energy_excited (float): Rigorous solution of the excited-state energy. + threshold (float): Regularization parameter (truncation lower bound). + return_vec (bool): If True, return value includes the corresponding eigenvector. + + Returns: + Union[float, Tuple[float, np.ndarray]]: The lowest eigenvalue (and its eigenvector if return_vec=True) + """ energies_excited = [] alpha_list = [] for i in range(1, len(self.qse_ops)): try: ret = self._solve_nlev_regularized_gen_eig( - h_eff, - s_mtrc, + h, + s, n_lev=i, threshold=threshold, return_vec=return_vec, @@ -217,7 +294,18 @@ def _solve_regularized_gen_eig_with_best_n_lev(self, h_eff, s_mtrc, true_energy_ else: return energies_excited[best_idx] - def execute_statistics(self, molecule): + def execute_statistics(self, molecule: MoleculeInfo) -> Tuple[float, float]: + """Evaluate mean absolute error and standard deviation of QSE calculation. + + Args: + molecule (MoleculeInfo): Target molecule + + Raises: + ValueError: Raised when the number of shots is too low and n_lev is too large. + + Returns: + Tuple[float, float]: Mean absolute error and standard deviation + """ self.h_ij, self.s_ij = self._get_h_s_ij(molecule.hamiltonian) h_eff_list = [] s_mtrc_list = []