diff --git a/src/simulators/matrix_product_state/matrix_product_state_internal.cpp b/src/simulators/matrix_product_state/matrix_product_state_internal.cpp index 779855f641..d512ef7e81 100644 --- a/src/simulators/matrix_product_state/matrix_product_state_internal.cpp +++ b/src/simulators/matrix_product_state/matrix_product_state_internal.cpp @@ -912,48 +912,71 @@ double MPS::expectation_value_internal(const reg_t &qubits, complex_t MPS::expectation_value_pauli(const reg_t &qubits, const std::string &matrices) const { reg_t internal_qubits = get_internal_qubits(qubits); - return expectation_value_pauli_internal(internal_qubits, matrices); -} -complex_t MPS::expectation_value_pauli_internal(const reg_t &qubits, const std::string &matrices) const { - reg_t sorted_qubits = qubits; - reg_t centralized_qubits = qubits; + // instead of computing the expectation value on the specified qubits, + // we find the min and max of these qubits, and compute the expectation value + // on all the qubits in between, inserting I matrices for those qubits + // that were not in the original vector "qubits". + // This enhancement was done for performance reasons + reg_t extended_qubits = internal_qubits; + + const auto min = std::min_element(begin(internal_qubits), end(internal_qubits)); + const auto max = std::max_element(begin(internal_qubits), end(internal_qubits)); + uint_t min_qubit = *min; + uint_t max_qubit = *max; + + // The number of qubits added to extended_qubits + uint_t num_Is = 0; + + // Add all the additional qubits at the end of the vector of extended_qubits + // The I matrices are added in expectation_value_pauli_internal, after they are reversed + for (uint_t i=min_qubit; i<=max_qubit; i++) { + auto itr = std::find(internal_qubits.begin(), internal_qubits.end(), i); + if (itr == internal_qubits.end()) { + extended_qubits.push_back(i); + num_Is++; + } + } + + return expectation_value_pauli_internal(extended_qubits, matrices, min_qubit, max_qubit, num_Is); +} - // if the qubits are not ordered, we can sort them, because the order doesn't matter +complex_t MPS::expectation_value_pauli_internal(const reg_t &qubits, + const std::string &matrices, + uint_t first_index, uint_t last_index, + uint_t num_Is) const { // when computing the expectation value. We only have to sort the pauli matrices // to be in the same ordering as the qubits - MPS temp_MPS; - MPS_with_new_indices(qubits, sorted_qubits, centralized_qubits, temp_MPS); - uint_t first_index = centralized_qubits.front(); - uint_t last_index = centralized_qubits.back(); - // Preliminary step - reverse the order of the matrices because + // Preliminary step - reverse the order of the matrices because // they are ordered in reverse to that of the qubits (in the interface) std::string reversed_matrices = matrices; reverse(reversed_matrices.begin(), reversed_matrices.end()); + for (uint_t i=0; i 0) { - left_tensor.mul_Gamma_by_left_Lambda(temp_MPS.lambda_reg_[first_index-1]); + left_tensor.mul_Gamma_by_left_Lambda(lambda_reg_[first_index-1]); } // The last gamma must be multiplied also by its right lambda. - // Here we handle the special case that we are calculating exp val + // Here we handle the special case that we are calculating exp val // on a single qubit // we need to mul every gamma by its right lambda if (first_index==last_index && first_index < num_qubits_-1) { - left_tensor.mul_Gamma_by_right_Lambda(temp_MPS.lambda_reg_[first_index]); + left_tensor.mul_Gamma_by_right_Lambda(lambda_reg_[first_index]); } - // Step 2 - prepare the dagger of left_tensor - MPS_Tensor left_tensor_dagger(AER::Utils::dagger(left_tensor.get_data(0)), + MPS_Tensor left_tensor_dagger(AER::Utils::dagger(left_tensor.get_data(0)), AER::Utils::dagger(left_tensor.get_data(1))); // Step 3 - Apply the gate to q0 left_tensor.apply_pauli(gate); @@ -966,38 +989,38 @@ complex_t MPS::expectation_value_pauli_internal(const reg_t &qubits, const std:: final_contract); for (uint_t qubit_num=first_index+1; qubit_num<=last_index; qubit_num++) { // Step 5 - multiply next Gamma by its left lambda (same as Step 1) - // next gamma has dimensions a0 x a1 x i - MPS_Tensor next_gamma = temp_MPS.q_reg_[qubit_num]; - next_gamma.mul_Gamma_by_left_Lambda(temp_MPS.lambda_reg_[qubit_num-1]); + // next gamma has dimensions a0 x a1 x i + MPS_Tensor next_gamma = q_reg_[qubit_num]; + next_gamma.mul_Gamma_by_left_Lambda(lambda_reg_[qubit_num-1]); // Last qubit must be multiplied by rightmost lambda if (qubit_num==last_index && qubit_num < num_qubits_-1) - next_gamma.mul_Gamma_by_right_Lambda(temp_MPS.lambda_reg_[qubit_num]); + next_gamma.mul_Gamma_by_right_Lambda(lambda_reg_[qubit_num]); // Step 6 - prepare the dagger of the next gamma (same as Step 2) // next_gamma_dagger has dimensions a1' x a0' x i - MPS_Tensor next_gamma_dagger(AER::Utils::dagger(next_gamma.get_data(0)), + MPS_Tensor next_gamma_dagger(AER::Utils::dagger(next_gamma.get_data(0)), AER::Utils::dagger(next_gamma.get_data(1))); - + // Step 7 - apply gate (same as Step 3) gate = sorted_matrices[qubit_num - first_index]; next_gamma.apply_pauli(gate); - + // Step 8 - contract final_contract from previous stage with next gamma over a1 // final_contract has dimensions a1 x a1, Gamma1 has dimensions a1 x a2 x i (where i=2) // result is a tensor of size a1 x a2 x i - MPS_Tensor next_contract(final_contract * next_gamma.get_data(0), + MPS_Tensor next_contract(final_contract * next_gamma.get_data(0), final_contract * next_gamma.get_data(1)); - // Step 9 - contract next_contract (a1 x a2 x i) + // Step 9 - contract next_contract (a1 x a2 x i) // with next_gamma_dagger (i x a2 x a1) (same as Step 4) // here we need to contract across two dimensions: a1 and i // result is a matrix of size a2 x a2 MPS_Tensor::contract_2_dimensions(next_gamma_dagger, next_contract, omp_threads_, - final_contract); + final_contract); } - + // Step 10 - contract over final matrix of size aN x aN // We need to contract the final matrix with itself // Compute this by taking the trace of final_contract diff --git a/src/simulators/matrix_product_state/matrix_product_state_internal.hpp b/src/simulators/matrix_product_state/matrix_product_state_internal.hpp index 2d4ff9d73e..91cae2db77 100644 --- a/src/simulators/matrix_product_state/matrix_product_state_internal.hpp +++ b/src/simulators/matrix_product_state/matrix_product_state_internal.hpp @@ -272,7 +272,9 @@ class MPS{ rvector_t trace_of_density_matrix(const reg_t &qubits) const; double expectation_value_internal(const reg_t &qubits, const cmatrix_t &M) const; - complex_t expectation_value_pauli_internal(const reg_t &qubits, const std::string &matrices) const; + complex_t expectation_value_pauli_internal(const reg_t &qubits, const std::string &matrices, + uint_t first_index, uint_t last_index, + uint_t num_Is) const; //---------------------------------------------------------------- // Function name: get_matrices_sizes