From e087712c5b5f60b8d97452f2006aa3bb17e32778 Mon Sep 17 00:00:00 2001 From: Ryan Babbush Date: Thu, 28 Sep 2017 19:12:19 -0700 Subject: [PATCH 1/4] Adding rdm measurement code --- src/openfermion/measurement/__init__.py | 16 ++ .../_equality_constraint_projection.py | 263 ++++++++++++++++++ .../_equality_constraint_projection_test.py | 150 ++++++++++ .../measurement/_rdm_equality_constraints.py | 115 ++++++++ .../_rdm_equality_constraints_test.py | 78 ++++++ 5 files changed, 622 insertions(+) create mode 100644 src/openfermion/measurement/__init__.py create mode 100644 src/openfermion/measurement/_equality_constraint_projection.py create mode 100644 src/openfermion/measurement/_equality_constraint_projection_test.py create mode 100644 src/openfermion/measurement/_rdm_equality_constraints.py create mode 100644 src/openfermion/measurement/_rdm_equality_constraints_test.py diff --git a/src/openfermion/measurement/__init__.py b/src/openfermion/measurement/__init__.py new file mode 100644 index 000000000..7ede95d70 --- /dev/null +++ b/src/openfermion/measurement/__init__.py @@ -0,0 +1,16 @@ +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +from ._equality_constraint_projection import apply_constraints + +from ._rdm_equality_constraints import (one_body_fermion_constraints, + two_body_fermion_constraints) diff --git a/src/openfermion/measurement/_equality_constraint_projection.py b/src/openfermion/measurement/_equality_constraint_projection.py new file mode 100644 index 000000000..315c912aa --- /dev/null +++ b/src/openfermion/measurement/_equality_constraint_projection.py @@ -0,0 +1,263 @@ +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +"""Module to reduce operator variance using equality RDM constraints.""" +import numpy +import scipy + +from openfermion.ops import FermionOperator, hermitian_conjugated +from openfermion.utils import count_qubits + +from _rdm_equality_constraints import two_body_fermion_constraints + + +def linearize_term(term, n_orbitals): + """Function to return integer index of term indices. + + Args: + term(tuple): The term indices of a one- or two-body FermionOperator. + n_orbitals(int): The number of orbitals in the simulation. + + Returns: + index(int): The index of the term. + """ + # Handle identity term. + if term == (): + return 0 + elif len(term) == 2: + # Handle one-body terms. + assert term[0][1] == 1 + assert term[1][1] == 0 + p = term[0][0] + q = term[1][0] + return 1 + p + q * n_orbitals + elif len(term) == 4: + # Handle two-body terms. + assert term[0][1] == 1 + assert term[1][1] == 1 + assert term[2][1] == 0 + assert term[3][1] == 0 + p = term[0][0] + q = term[1][0] + r = term[2][0] + s = term[3][0] + return (1 + n_orbitals ** 2 + + p + + q * n_orbitals + + r * n_orbitals ** 2 + + s * n_orbitals ** 3) + + +def unlinearize_term(index, n_orbitals): + """Function to return integer index of term indices. + + Args: + index(int): The index of the term. + n_orbitals(int): The number of orbitals in the simulation. + + Returns: + term(tuple): The term indices of a one- or two-body FermionOperator. + """ + # Handle identity term. + if not index: + return (()) + elif (0 < index < 1 + n_orbitals ** 2): + # Handle one-body terms. + shift = 1 + new_index = index - shift + q = new_index // n_orbitals + p = new_index - q * n_orbitals + assert index == shift + p + q * n_orbitals + return ((p, 1), (q, 0)) + else: + # Handle two-body terms. + shift = 1 + n_orbitals ** 2 + new_index = index - shift + s = new_index // n_orbitals ** 3 + r = (new_index - s * n_orbitals ** 3) // n_orbitals ** 2 + q = (new_index - s * n_orbitals ** 3 - + r * n_orbitals ** 2) // n_orbitals + p = (new_index - q * n_orbitals - + r * n_orbitals ** 2 - s * n_orbitals ** 3) + assert index == (shift + p + q * n_orbitals + + r * n_orbitals ** 2 + s * n_orbitals ** 3) + return ((p, 1), (q, 1), (r, 0), (s, 0)) + + +def constraint_matrix(n_orbitals, n_fermions): + """Function to generate matrix of constraints. + + Args: + n_orbitals(int): The number of orbitals in the simulation. + n_fermions(int): The number of particles in the simulation. + + Returns: + constraint_matrix(scipy.sparse.coo_matrix): The matrix of constraints. + """ + # Very inefficiently count constraints. + n_constraints = 0 + for constraint in two_body_fermion_constraints(n_orbitals, n_fermions): + n_constraints += 1 + + # Initialize constraint matrix. + n_terms = 1 + n_orbitals ** 2 + n_orbitals ** 4 + constraint_matrix = scipy.sparse.dok_matrix((n_constraints, n_terms)) + + # Populate constraint matrix. + constraint_number = 0 + for constraint in two_body_fermion_constraints(n_orbitals, n_fermions): + for term, coefficient in constraint.terms.items(): + term_index = linearize_term(term, n_orbitals) + constraint_matrix[constraint_number, term_index] = coefficient + constraint_number += 1 + return constraint_matrix + + +def operator_to_vector(operator): + """Function to map operator to vector. + + Args: + operator(FermionOperator): FermionOperator with only 1- and 2-body + terms that we wish to vectorize. + + Returns: + vectorized_operator(numpy.array): Vector of term coefficients. + """ + n_orbitals = count_qubits(operator) + n_terms = 1 + n_orbitals ** 2 + n_orbitals ** 4 + vectorized_operator = numpy.zeros(n_terms, float) + for term, coefficient in operator.terms.items(): + term_index = linearize_term(term, n_orbitals) + vectorized_operator[term_index] = coefficient + return vectorized_operator + + +def vector_to_operator(vector, n_orbitals): + """Function to map vector to operator. + + Args: + vectorized_operator(numpy.array): Vector of term coefficients. + + Returns: + operator(FermionOperator): FermionOperator with only 1- and 2-body + terms that we wish to vectorize. + """ + operator = FermionOperator() + for index, coefficient in enumerate(vector): + term = unlinearize_term(index, n_orbitals) + operator += FermionOperator(term, coefficient) + return operator + + +def apply_constraints(operator, n_fermions, use_scipy=True): + """Function to use linear programming to apply constraints. + + Args: + operator(FermionOperator): FermionOperator with only 1- and 2-body + terms that we wish to vectorize. + n_fermions(int): The number of particles in the simulation. + use_scipy(bool): Whether to use scipy (True) or cvxopt (False). + + Returns: + modified_operator(FermionOperator): The operator with reduced norm + that has been modified with equality constraints. + """ + # Get constraint matrix. + n_orbitals = count_qubits(operator) + constraints = constraint_matrix(n_orbitals, n_fermions) + n_constraints, n_terms = constraints.get_shape() + + # Get vectorized operator. + vectorized_operator = operator_to_vector(operator) + initial_bound = numpy.sum(numpy.absolute(vectorized_operator[1::])) ** 2 + print('Initial bound on measurements is %f.' % initial_bound) + + # Get linear programming coefficient vector. + n_variables = n_constraints + n_terms + lp_vector = numpy.zeros(n_variables, float) + lp_vector[-n_terms:] = 1. + + # Get linear programming constraint matrix. + lp_constraint_matrix = scipy.sparse.dok_matrix((2 * n_terms, n_variables)) + for (i, j), value in constraints.items(): + if j: + lp_constraint_matrix[j, i] = value + lp_constraint_matrix[n_terms + j, i] = -value + for i in range(n_terms): + lp_constraint_matrix[i, n_constraints + i] = -1. + lp_constraint_matrix[n_terms + i, n_constraints + i] = -1. + + # Get linear programming constraint vector. + lp_constraint_vector = numpy.zeros(2 * n_terms, float) + lp_constraint_vector[:n_terms] = vectorized_operator + lp_constraint_vector[-n_terms:] = -vectorized_operator + + # Perform linear programming. + print('Starting linear programming.') + if use_scipy: + options = {'maxiter': int(1e6)} + bound = n_constraints * [(None, None)] + n_terms * [(0, None)] + solution = scipy.optimize.linprog(c=lp_vector, + A_ub=lp_constraint_matrix.toarray(), + b_ub=lp_constraint_vector, + bounds=bound, + options=options) + + # Analyze results. + print(solution['message']) + assert solution['success'] + solution_vector = solution['x'] + objective = solution['fun'] ** 2 + print('Program terminated after %i iterations.' % solution['nit']) + + else: + # Convert to CVXOpt sparse matrix. + from cvxopt import matrix, solvers, spmatrix + lp_vector = matrix(lp_vector) + lp_constraint_matrix = lp_constraint_matrix.tocoo() + lp_constraint_matrix = spmatrix(lp_constraint_matrix.data, + lp_constraint_matrix.row.tolist(), + lp_constraint_matrix.col.tolist()) + lp_constraint_vector = matrix(lp_constraint_vector) + + # Run linear programming. + solution = solvers.lp(c=lp_vector, + G=lp_constraint_matrix, + h=lp_constraint_vector, + solver='glpk') + + # Analyze results. + print(solution['status']) + solution_vector = numpy.array(solution['x']).transpose()[0] + + # Alternative bound. + residuals = solution_vector[-n_terms:] + alternative_bound = numpy.sum(numpy.absolute(residuals[1::])) ** 2 + print('Bound implied by solution vector is %f.' % alternative_bound) + + # Make sure residuals are positive. + for residual in residuals: + assert residual > -1e-6 + + # Get bound on updated Hamiltonian. + weights = solution_vector[:n_constraints] + final_vectorized_operator = (vectorized_operator - + constraints.transpose() * weights) + final_bound = numpy.sum( + numpy.absolute(final_vectorized_operator[1::])) ** 2 + print('Actual bound determined is %f.' % final_bound) + + # Return modified operator. + modified_operator = vector_to_operator( + final_vectorized_operator, n_orbitals) + return (modified_operator + + hermitian_conjugated(modified_operator)) / 2. diff --git a/src/openfermion/measurement/_equality_constraint_projection_test.py b/src/openfermion/measurement/_equality_constraint_projection_test.py new file mode 100644 index 000000000..419b18718 --- /dev/null +++ b/src/openfermion/measurement/_equality_constraint_projection_test.py @@ -0,0 +1,150 @@ +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +"""Tests for _variance_reduction.py""" +import numpy +import os +import unittest + +from _equality_constraint_projection import (apply_constraints, + constraint_matrix, + linearize_term, + operator_to_vector, + unlinearize_term, + vector_to_operator) + +from openfermion.config import THIS_DIRECTORY +from openfermion.transforms import get_fermion_operator, get_sparse_operator +from openfermion.utils import (count_qubits, expectation, + MolecularData, get_ground_state) +from openfermion.utils._sparse_tools import (jw_number_restrict_operator, + sparse_eigenspectrum) + + +class EqualityConstraintProjectionTest(unittest.TestCase): + + def setUp(self): + + # Set up molecule. + n_atoms = 2 + geometry = [('H', (0., 0., 0.)), ('H', (0., 0., 0.7414))] + basis = 'sto-3g' + multiplicity = 1 + filename = os.path.join(THIS_DIRECTORY, 'data', + 'H2_sto-3g_singlet_0.7414') + molecule = MolecularData( + geometry, basis, multiplicity, filename=filename) + molecule.load() + self.n_fermions = molecule.n_electrons + self.n_orbitals = molecule.n_qubits + + # Get molecular Hamiltonian. + molecular_hamiltonian = molecule.get_molecular_hamiltonian() + self.fermion_hamiltonian = get_fermion_operator(molecular_hamiltonian) + + def test_linearize_term(self): + past_terms = set() + for term, coefficient in self.fermion_hamiltonian.terms.items(): + index = linearize_term(term, self.n_orbitals) + self.assertTrue(isinstance(index, int)) + self.assertFalse(index in past_terms) + past_terms.add(index) + + def test_unlinearize_term_consistency(self): + for term, coefficient in self.fermion_hamiltonian.terms.items(): + index = linearize_term(term, self.n_orbitals) + new_term = unlinearize_term(index, self.n_orbitals) + self.assertEqual(term, new_term) + + def test_operator_to_vector_consistency(self): + vector = operator_to_vector(self.fermion_hamiltonian) + operator = vector_to_operator(vector, self.n_orbitals) + magnitude = 0. + difference = operator - self.fermion_hamiltonian + for term, coefficient in difference.terms.items(): + magnitude += abs(coefficient) + self.assertAlmostEqual(0, magnitude) + + def test_constraint_matrix(self): + + # Randomly project operator with constraints. + numpy.random.seed(8) + constraints = constraint_matrix(self.n_orbitals, self.n_fermions) + n_constraints, n_terms = constraints.shape + self.assertEqual( + 1 + self.n_orbitals ** 2 + self.n_orbitals ** 4, n_terms) + random_weights = numpy.random.randn(n_constraints) + vectorized_operator = operator_to_vector(self.fermion_hamiltonian) + modification_vector = constraints.transpose() * random_weights + new_operator_vector = vectorized_operator + modification_vector + modified_operator = vector_to_operator( + new_operator_vector, self.n_orbitals) + + # Map both to sparse matrix under Jordan-Wigner. + sparse_original = get_sparse_operator(self.fermion_hamiltonian) + sparse_modified = get_sparse_operator(modified_operator) + + # Check expectation value. + energy, wavefunction = get_ground_state(sparse_original) + modified_energy = expectation(sparse_modified, wavefunction) + self.assertAlmostEqual(modified_energy, energy) + + def test_apply_constraints(self): + + # Get norm of original operator. + original_norm = 0. + for term, coefficient in self.fermion_hamiltonian.terms.items(): + if term != (): + original_norm += abs(coefficient) + + # Get modified operator. + modified_operator = apply_constraints( + self.fermion_hamiltonian, self.n_fermions) + modified_operator.compress() + + # Get norm of modified operator. + modified_norm = 0. + for term, coefficient in modified_operator.terms.items(): + if term != (): + modified_norm += abs(coefficient) + self.assertTrue(modified_norm < original_norm) + + # Map both to sparse matrix under Jordan-Wigner. + sparse_original = get_sparse_operator(self.fermion_hamiltonian) + sparse_modified = get_sparse_operator(modified_operator) + + # Check spectra. + sparse_original = jw_number_restrict_operator(sparse_original, + self.n_fermions) + sparse_modified = jw_number_restrict_operator(sparse_modified, + self.n_fermions) + original_spectrum = sparse_eigenspectrum(sparse_original) + modified_spectrum = sparse_eigenspectrum(sparse_modified) + spectral_deviation = numpy.amax(numpy.absolute( + original_spectrum - modified_spectrum)) + self.assertAlmostEqual(spectral_deviation, 0.) + + # Check expectation value. + energy, wavefunction = get_ground_state(sparse_original) + modified_energy = expectation(sparse_modified, wavefunction) + self.assertAlmostEqual(modified_energy, energy) + + # Test consistency with cvxopt. + scipy_operator = apply_constraints( + self.fermion_hamiltonian, self.n_fermions, use_scipy=False) + + # Get norm. + scipy_norm = 0. + for term, coefficient in scipy_operator.terms.items(): + if term != (): + scipy_norm += abs(coefficient) + self.assertAlmostEqual(scipy_norm, modified_norm) diff --git a/src/openfermion/measurement/_rdm_equality_constraints.py b/src/openfermion/measurement/_rdm_equality_constraints.py new file mode 100644 index 000000000..2291aa188 --- /dev/null +++ b/src/openfermion/measurement/_rdm_equality_constraints.py @@ -0,0 +1,115 @@ +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +"""Constraints on fermionic reduced density matrices""" +from openfermion.ops import FermionOperator + + +def one_body_fermion_constraints(n_orbitals, n_fermions): + """Generates one-body positivity constraints on fermionic RDMs. + + The specific constraints implemented are known positivity constraints + on the one-fermion reduced density matrices. Constraints are generated + in the form of FermionOperators whose expectation value is known to be + zero for any N-Representable state. Generators are used for efficiency. + + Args: + n_orbitals(int): number of spin-orbitals on which operators act. + n_fermions(int): number of fermions in the system. + + Yields: + Constraint is a FermionOperator with zero expectation value. + """ + # One-RDM trace condition. + constraint_operator = FermionOperator() + for i in range(n_orbitals): + constraint_operator += FermionOperator(((i, 1), (i, 0))) + if len(constraint_operator.terms): + constraint_operator -= FermionOperator((), n_fermions) + yield constraint_operator + + # One-RDM Hermiticity condition. + for i in range(n_orbitals): + for j in range(i + 1, n_orbitals): + constraint_operator = FermionOperator(((i, 1), (j, 0))) + constraint_operator -= FermionOperator(((j, 1), (i, 0))) + if len(constraint_operator.terms): + yield constraint_operator + + +def two_body_fermion_constraints(n_orbitals, n_fermions): + """Generates two-body positivity constraints on fermionic RDMs. + + The specific constraints implemented are known positivity constraints + on the two-fermion reduced density matrices. Constraints are generated + in the form of FermionOperators whose expectation value is known to be + zero for any N-Representable state. Generators are used for efficiency. + + Args: + n_orbitals(int): number of spin-orbitals on which operators act. + n_fermions(int): number of fermions in the system. + + Yields: + Constraint is a FermionOperator with zero expectation value. + """ + # Two-body trace condition. + constraint_operator = FermionOperator() + for i in range(n_orbitals): + for j in range(n_orbitals): + constraint_operator += FermionOperator( + ((i, 1), (j, 1), (j, 0), (i, 0))) + if len(constraint_operator.terms): + constraint_operator -= FermionOperator( + (), n_fermions * (n_fermions - 1)) + yield constraint_operator + + # Two-body Hermiticity condition. + for ij in range(n_orbitals ** 2): + i, j = (ij / n_orbitals), (ij % n_orbitals) + for kl in range(ij + 1, n_orbitals ** 2): + k, l = (kl / n_orbitals), (kl % n_orbitals) + constraint_operator = FermionOperator( + ((i, 1), (j, 1), (l, 0), (k, 0))) + constraint_operator -= FermionOperator( + ((k, 1), (l, 1), (j, 0), (i, 0))) + if len(constraint_operator.terms): + yield constraint_operator + + # Contraction to One-RDM from Two-RDM. + for i in range(n_orbitals): + for j in range(n_orbitals): + constraint_operator = FermionOperator() + for p in range(n_orbitals): + constraint_operator += FermionOperator( + ((i, 1), (p, 1), (p, 0), (j, 0))) + constraint_operator += FermionOperator( + ((i, 1), (j, 0)), -(n_fermions - 1)) + if len(constraint_operator.terms): + yield constraint_operator + + # Linear relations between two-particle matrices. + for ij in range(n_orbitals ** 2): + i, j = (ij / n_orbitals), (ij % n_orbitals) + for kl in range(ij, n_orbitals ** 2): + k, l = (kl / n_orbitals), (kl % n_orbitals) + + # G-matrix condition. + constraint_operator = FermionOperator( + ((i, 1), (k, 0)), 1.0 * (j == l)) + constraint_operator += FermionOperator( + ((i, 1), (l, 1), (k, 0), (j, 0))) + constraint_operator += FermionOperator( + ((i, 1), (l, 1), (j, 0), (k, 0))) + constraint_operator -= FermionOperator( + ((i, 1), (k, 0)), 1.0 * (j == l)) + if len(constraint_operator.terms): + yield constraint_operator diff --git a/src/openfermion/measurement/_rdm_equality_constraints_test.py b/src/openfermion/measurement/_rdm_equality_constraints_test.py new file mode 100644 index 000000000..c7e2f5328 --- /dev/null +++ b/src/openfermion/measurement/_rdm_equality_constraints_test.py @@ -0,0 +1,78 @@ +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +"""Tests for _rdm_equality_constraints.py""" +import unittest +import os + +from _rdm_equality_constraints import (one_body_fermion_constraints, + two_body_fermion_constraints) + +from openfermion.config import THIS_DIRECTORY +from openfermion.ops import FermionOperator, hermitian_conjugated +from openfermion.transforms import get_interaction_operator +from openfermion.utils import MolecularData + + +class FermionConstraintsTest(unittest.TestCase): + + def setUp(self): + + # Setup. + n_atoms = 2 + geometry = [('H', (0., 0., 0.)), ('H', (0., 0., 0.7414))] + basis = 'sto-3g' + multiplicity = 1 + filename = os.path.join(THIS_DIRECTORY, 'data', + 'H2_sto-3g_singlet_0.7414') + molecule = MolecularData( + geometry, basis, multiplicity, filename=filename) + molecule.load() + self.n_fermions = molecule.n_electrons + self.n_orbitals = molecule.n_qubits + + # Get molecular Hamiltonian. + self.molecular_hamiltonian = molecule.get_molecular_hamiltonian() + self.fci_rdm = molecule.get_molecular_rdm(use_fci=1) + + def test_one_body_constraints(self): + for constraint in one_body_fermion_constraints( + self.n_orbitals, self.n_fermions): + interaction_operator = get_interaction_operator( + constraint, self.n_orbitals) + constraint_value = self.fci_rdm.expectation(interaction_operator) + self.assertAlmostEqual(constraint_value, 0.) + for term, coefficient in constraint.terms.items(): + if len(term) == 2: + self.assertTrue(term[0][1]) + self.assertFalse(term[1][1]) + else: + self.assertEqual(term, ()) + + def test_two_body_constraints(self): + for constraint in two_body_fermion_constraints( + self.n_orbitals, self.n_fermions): + interaction_operator = get_interaction_operator( + constraint, self.n_orbitals) + constraint_value = self.fci_rdm.expectation(interaction_operator) + self.assertAlmostEqual(constraint_value, 0.) + for term, coefficient in constraint.terms.items(): + if len(term) == 2: + self.assertTrue(term[0][1]) + self.assertFalse(term[1][1]) + elif len(term) == 4: + self.assertTrue(term[0][1]) + self.assertTrue(term[1][1]) + self.assertFalse(term[2][1]) + self.assertFalse(term[3][1]) + else: + self.assertEqual(term, ()) From 0bb3a1e8afdc57a844621fe98612e18c5f1aee28 Mon Sep 17 00:00:00 2001 From: Ryan Babbush Date: Thu, 28 Sep 2017 19:27:01 -0700 Subject: [PATCH 2/4] changing pull request to only include constraint code for now --- src/openfermion/measurement/__init__.py | 16 -- .../_equality_constraint_projection.py | 263 ------------------ .../_equality_constraint_projection_test.py | 150 ---------- src/openfermion/utils/__init__.py | 4 + .../_rdm_equality_constraints.py | 0 .../_rdm_equality_constraints_test.py | 8 +- src/openfermion/utils/_sparse_tools.py | 2 +- 7 files changed, 8 insertions(+), 435 deletions(-) delete mode 100644 src/openfermion/measurement/__init__.py delete mode 100644 src/openfermion/measurement/_equality_constraint_projection.py delete mode 100644 src/openfermion/measurement/_equality_constraint_projection_test.py rename src/openfermion/{measurement => utils}/_rdm_equality_constraints.py (100%) rename src/openfermion/{measurement => utils}/_rdm_equality_constraints_test.py (92%) diff --git a/src/openfermion/measurement/__init__.py b/src/openfermion/measurement/__init__.py deleted file mode 100644 index 7ede95d70..000000000 --- a/src/openfermion/measurement/__init__.py +++ /dev/null @@ -1,16 +0,0 @@ -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. - -from ._equality_constraint_projection import apply_constraints - -from ._rdm_equality_constraints import (one_body_fermion_constraints, - two_body_fermion_constraints) diff --git a/src/openfermion/measurement/_equality_constraint_projection.py b/src/openfermion/measurement/_equality_constraint_projection.py deleted file mode 100644 index 315c912aa..000000000 --- a/src/openfermion/measurement/_equality_constraint_projection.py +++ /dev/null @@ -1,263 +0,0 @@ -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. - -"""Module to reduce operator variance using equality RDM constraints.""" -import numpy -import scipy - -from openfermion.ops import FermionOperator, hermitian_conjugated -from openfermion.utils import count_qubits - -from _rdm_equality_constraints import two_body_fermion_constraints - - -def linearize_term(term, n_orbitals): - """Function to return integer index of term indices. - - Args: - term(tuple): The term indices of a one- or two-body FermionOperator. - n_orbitals(int): The number of orbitals in the simulation. - - Returns: - index(int): The index of the term. - """ - # Handle identity term. - if term == (): - return 0 - elif len(term) == 2: - # Handle one-body terms. - assert term[0][1] == 1 - assert term[1][1] == 0 - p = term[0][0] - q = term[1][0] - return 1 + p + q * n_orbitals - elif len(term) == 4: - # Handle two-body terms. - assert term[0][1] == 1 - assert term[1][1] == 1 - assert term[2][1] == 0 - assert term[3][1] == 0 - p = term[0][0] - q = term[1][0] - r = term[2][0] - s = term[3][0] - return (1 + n_orbitals ** 2 + - p + - q * n_orbitals + - r * n_orbitals ** 2 + - s * n_orbitals ** 3) - - -def unlinearize_term(index, n_orbitals): - """Function to return integer index of term indices. - - Args: - index(int): The index of the term. - n_orbitals(int): The number of orbitals in the simulation. - - Returns: - term(tuple): The term indices of a one- or two-body FermionOperator. - """ - # Handle identity term. - if not index: - return (()) - elif (0 < index < 1 + n_orbitals ** 2): - # Handle one-body terms. - shift = 1 - new_index = index - shift - q = new_index // n_orbitals - p = new_index - q * n_orbitals - assert index == shift + p + q * n_orbitals - return ((p, 1), (q, 0)) - else: - # Handle two-body terms. - shift = 1 + n_orbitals ** 2 - new_index = index - shift - s = new_index // n_orbitals ** 3 - r = (new_index - s * n_orbitals ** 3) // n_orbitals ** 2 - q = (new_index - s * n_orbitals ** 3 - - r * n_orbitals ** 2) // n_orbitals - p = (new_index - q * n_orbitals - - r * n_orbitals ** 2 - s * n_orbitals ** 3) - assert index == (shift + p + q * n_orbitals + - r * n_orbitals ** 2 + s * n_orbitals ** 3) - return ((p, 1), (q, 1), (r, 0), (s, 0)) - - -def constraint_matrix(n_orbitals, n_fermions): - """Function to generate matrix of constraints. - - Args: - n_orbitals(int): The number of orbitals in the simulation. - n_fermions(int): The number of particles in the simulation. - - Returns: - constraint_matrix(scipy.sparse.coo_matrix): The matrix of constraints. - """ - # Very inefficiently count constraints. - n_constraints = 0 - for constraint in two_body_fermion_constraints(n_orbitals, n_fermions): - n_constraints += 1 - - # Initialize constraint matrix. - n_terms = 1 + n_orbitals ** 2 + n_orbitals ** 4 - constraint_matrix = scipy.sparse.dok_matrix((n_constraints, n_terms)) - - # Populate constraint matrix. - constraint_number = 0 - for constraint in two_body_fermion_constraints(n_orbitals, n_fermions): - for term, coefficient in constraint.terms.items(): - term_index = linearize_term(term, n_orbitals) - constraint_matrix[constraint_number, term_index] = coefficient - constraint_number += 1 - return constraint_matrix - - -def operator_to_vector(operator): - """Function to map operator to vector. - - Args: - operator(FermionOperator): FermionOperator with only 1- and 2-body - terms that we wish to vectorize. - - Returns: - vectorized_operator(numpy.array): Vector of term coefficients. - """ - n_orbitals = count_qubits(operator) - n_terms = 1 + n_orbitals ** 2 + n_orbitals ** 4 - vectorized_operator = numpy.zeros(n_terms, float) - for term, coefficient in operator.terms.items(): - term_index = linearize_term(term, n_orbitals) - vectorized_operator[term_index] = coefficient - return vectorized_operator - - -def vector_to_operator(vector, n_orbitals): - """Function to map vector to operator. - - Args: - vectorized_operator(numpy.array): Vector of term coefficients. - - Returns: - operator(FermionOperator): FermionOperator with only 1- and 2-body - terms that we wish to vectorize. - """ - operator = FermionOperator() - for index, coefficient in enumerate(vector): - term = unlinearize_term(index, n_orbitals) - operator += FermionOperator(term, coefficient) - return operator - - -def apply_constraints(operator, n_fermions, use_scipy=True): - """Function to use linear programming to apply constraints. - - Args: - operator(FermionOperator): FermionOperator with only 1- and 2-body - terms that we wish to vectorize. - n_fermions(int): The number of particles in the simulation. - use_scipy(bool): Whether to use scipy (True) or cvxopt (False). - - Returns: - modified_operator(FermionOperator): The operator with reduced norm - that has been modified with equality constraints. - """ - # Get constraint matrix. - n_orbitals = count_qubits(operator) - constraints = constraint_matrix(n_orbitals, n_fermions) - n_constraints, n_terms = constraints.get_shape() - - # Get vectorized operator. - vectorized_operator = operator_to_vector(operator) - initial_bound = numpy.sum(numpy.absolute(vectorized_operator[1::])) ** 2 - print('Initial bound on measurements is %f.' % initial_bound) - - # Get linear programming coefficient vector. - n_variables = n_constraints + n_terms - lp_vector = numpy.zeros(n_variables, float) - lp_vector[-n_terms:] = 1. - - # Get linear programming constraint matrix. - lp_constraint_matrix = scipy.sparse.dok_matrix((2 * n_terms, n_variables)) - for (i, j), value in constraints.items(): - if j: - lp_constraint_matrix[j, i] = value - lp_constraint_matrix[n_terms + j, i] = -value - for i in range(n_terms): - lp_constraint_matrix[i, n_constraints + i] = -1. - lp_constraint_matrix[n_terms + i, n_constraints + i] = -1. - - # Get linear programming constraint vector. - lp_constraint_vector = numpy.zeros(2 * n_terms, float) - lp_constraint_vector[:n_terms] = vectorized_operator - lp_constraint_vector[-n_terms:] = -vectorized_operator - - # Perform linear programming. - print('Starting linear programming.') - if use_scipy: - options = {'maxiter': int(1e6)} - bound = n_constraints * [(None, None)] + n_terms * [(0, None)] - solution = scipy.optimize.linprog(c=lp_vector, - A_ub=lp_constraint_matrix.toarray(), - b_ub=lp_constraint_vector, - bounds=bound, - options=options) - - # Analyze results. - print(solution['message']) - assert solution['success'] - solution_vector = solution['x'] - objective = solution['fun'] ** 2 - print('Program terminated after %i iterations.' % solution['nit']) - - else: - # Convert to CVXOpt sparse matrix. - from cvxopt import matrix, solvers, spmatrix - lp_vector = matrix(lp_vector) - lp_constraint_matrix = lp_constraint_matrix.tocoo() - lp_constraint_matrix = spmatrix(lp_constraint_matrix.data, - lp_constraint_matrix.row.tolist(), - lp_constraint_matrix.col.tolist()) - lp_constraint_vector = matrix(lp_constraint_vector) - - # Run linear programming. - solution = solvers.lp(c=lp_vector, - G=lp_constraint_matrix, - h=lp_constraint_vector, - solver='glpk') - - # Analyze results. - print(solution['status']) - solution_vector = numpy.array(solution['x']).transpose()[0] - - # Alternative bound. - residuals = solution_vector[-n_terms:] - alternative_bound = numpy.sum(numpy.absolute(residuals[1::])) ** 2 - print('Bound implied by solution vector is %f.' % alternative_bound) - - # Make sure residuals are positive. - for residual in residuals: - assert residual > -1e-6 - - # Get bound on updated Hamiltonian. - weights = solution_vector[:n_constraints] - final_vectorized_operator = (vectorized_operator - - constraints.transpose() * weights) - final_bound = numpy.sum( - numpy.absolute(final_vectorized_operator[1::])) ** 2 - print('Actual bound determined is %f.' % final_bound) - - # Return modified operator. - modified_operator = vector_to_operator( - final_vectorized_operator, n_orbitals) - return (modified_operator + - hermitian_conjugated(modified_operator)) / 2. diff --git a/src/openfermion/measurement/_equality_constraint_projection_test.py b/src/openfermion/measurement/_equality_constraint_projection_test.py deleted file mode 100644 index 419b18718..000000000 --- a/src/openfermion/measurement/_equality_constraint_projection_test.py +++ /dev/null @@ -1,150 +0,0 @@ -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. - -"""Tests for _variance_reduction.py""" -import numpy -import os -import unittest - -from _equality_constraint_projection import (apply_constraints, - constraint_matrix, - linearize_term, - operator_to_vector, - unlinearize_term, - vector_to_operator) - -from openfermion.config import THIS_DIRECTORY -from openfermion.transforms import get_fermion_operator, get_sparse_operator -from openfermion.utils import (count_qubits, expectation, - MolecularData, get_ground_state) -from openfermion.utils._sparse_tools import (jw_number_restrict_operator, - sparse_eigenspectrum) - - -class EqualityConstraintProjectionTest(unittest.TestCase): - - def setUp(self): - - # Set up molecule. - n_atoms = 2 - geometry = [('H', (0., 0., 0.)), ('H', (0., 0., 0.7414))] - basis = 'sto-3g' - multiplicity = 1 - filename = os.path.join(THIS_DIRECTORY, 'data', - 'H2_sto-3g_singlet_0.7414') - molecule = MolecularData( - geometry, basis, multiplicity, filename=filename) - molecule.load() - self.n_fermions = molecule.n_electrons - self.n_orbitals = molecule.n_qubits - - # Get molecular Hamiltonian. - molecular_hamiltonian = molecule.get_molecular_hamiltonian() - self.fermion_hamiltonian = get_fermion_operator(molecular_hamiltonian) - - def test_linearize_term(self): - past_terms = set() - for term, coefficient in self.fermion_hamiltonian.terms.items(): - index = linearize_term(term, self.n_orbitals) - self.assertTrue(isinstance(index, int)) - self.assertFalse(index in past_terms) - past_terms.add(index) - - def test_unlinearize_term_consistency(self): - for term, coefficient in self.fermion_hamiltonian.terms.items(): - index = linearize_term(term, self.n_orbitals) - new_term = unlinearize_term(index, self.n_orbitals) - self.assertEqual(term, new_term) - - def test_operator_to_vector_consistency(self): - vector = operator_to_vector(self.fermion_hamiltonian) - operator = vector_to_operator(vector, self.n_orbitals) - magnitude = 0. - difference = operator - self.fermion_hamiltonian - for term, coefficient in difference.terms.items(): - magnitude += abs(coefficient) - self.assertAlmostEqual(0, magnitude) - - def test_constraint_matrix(self): - - # Randomly project operator with constraints. - numpy.random.seed(8) - constraints = constraint_matrix(self.n_orbitals, self.n_fermions) - n_constraints, n_terms = constraints.shape - self.assertEqual( - 1 + self.n_orbitals ** 2 + self.n_orbitals ** 4, n_terms) - random_weights = numpy.random.randn(n_constraints) - vectorized_operator = operator_to_vector(self.fermion_hamiltonian) - modification_vector = constraints.transpose() * random_weights - new_operator_vector = vectorized_operator + modification_vector - modified_operator = vector_to_operator( - new_operator_vector, self.n_orbitals) - - # Map both to sparse matrix under Jordan-Wigner. - sparse_original = get_sparse_operator(self.fermion_hamiltonian) - sparse_modified = get_sparse_operator(modified_operator) - - # Check expectation value. - energy, wavefunction = get_ground_state(sparse_original) - modified_energy = expectation(sparse_modified, wavefunction) - self.assertAlmostEqual(modified_energy, energy) - - def test_apply_constraints(self): - - # Get norm of original operator. - original_norm = 0. - for term, coefficient in self.fermion_hamiltonian.terms.items(): - if term != (): - original_norm += abs(coefficient) - - # Get modified operator. - modified_operator = apply_constraints( - self.fermion_hamiltonian, self.n_fermions) - modified_operator.compress() - - # Get norm of modified operator. - modified_norm = 0. - for term, coefficient in modified_operator.terms.items(): - if term != (): - modified_norm += abs(coefficient) - self.assertTrue(modified_norm < original_norm) - - # Map both to sparse matrix under Jordan-Wigner. - sparse_original = get_sparse_operator(self.fermion_hamiltonian) - sparse_modified = get_sparse_operator(modified_operator) - - # Check spectra. - sparse_original = jw_number_restrict_operator(sparse_original, - self.n_fermions) - sparse_modified = jw_number_restrict_operator(sparse_modified, - self.n_fermions) - original_spectrum = sparse_eigenspectrum(sparse_original) - modified_spectrum = sparse_eigenspectrum(sparse_modified) - spectral_deviation = numpy.amax(numpy.absolute( - original_spectrum - modified_spectrum)) - self.assertAlmostEqual(spectral_deviation, 0.) - - # Check expectation value. - energy, wavefunction = get_ground_state(sparse_original) - modified_energy = expectation(sparse_modified, wavefunction) - self.assertAlmostEqual(modified_energy, energy) - - # Test consistency with cvxopt. - scipy_operator = apply_constraints( - self.fermion_hamiltonian, self.n_fermions, use_scipy=False) - - # Get norm. - scipy_norm = 0. - for term, coefficient in scipy_operator.terms.items(): - if term != (): - scipy_norm += abs(coefficient) - self.assertAlmostEqual(scipy_norm, modified_norm) diff --git a/src/openfermion/utils/__init__.py b/src/openfermion/utils/__init__.py index 9c5610e30..335db0506 100644 --- a/src/openfermion/utils/__init__.py +++ b/src/openfermion/utils/__init__.py @@ -43,6 +43,9 @@ jordan_wigner_dual_basis_hamiltonian, wigner_seitz_length_scale) +from ._rdm_equality_constraints import (one_body_fermion_constraints, + two_body_fermion_constraints) + from ._sparse_tools import (expectation, expectation_computational_basis_state, get_density_matrix, @@ -51,6 +54,7 @@ is_hermitian, jordan_wigner_sparse, jw_hartree_fock_state, + jw_number_restrict_operator, qubit_operator_sparse, sparse_eigenspectrum) diff --git a/src/openfermion/measurement/_rdm_equality_constraints.py b/src/openfermion/utils/_rdm_equality_constraints.py similarity index 100% rename from src/openfermion/measurement/_rdm_equality_constraints.py rename to src/openfermion/utils/_rdm_equality_constraints.py diff --git a/src/openfermion/measurement/_rdm_equality_constraints_test.py b/src/openfermion/utils/_rdm_equality_constraints_test.py similarity index 92% rename from src/openfermion/measurement/_rdm_equality_constraints_test.py rename to src/openfermion/utils/_rdm_equality_constraints_test.py index c7e2f5328..d3f4c8642 100644 --- a/src/openfermion/measurement/_rdm_equality_constraints_test.py +++ b/src/openfermion/utils/_rdm_equality_constraints_test.py @@ -14,13 +14,11 @@ import unittest import os -from _rdm_equality_constraints import (one_body_fermion_constraints, - two_body_fermion_constraints) - from openfermion.config import THIS_DIRECTORY -from openfermion.ops import FermionOperator, hermitian_conjugated from openfermion.transforms import get_interaction_operator -from openfermion.utils import MolecularData +from openfermion.utils import (MolecularData, + one_body_fermion_constraints, + two_body_fermion_constraints) class FermionConstraintsTest(unittest.TestCase): diff --git a/src/openfermion/utils/_sparse_tools.py b/src/openfermion/utils/_sparse_tools.py index c7f78e7ff..b50927b47 100644 --- a/src/openfermion/utils/_sparse_tools.py +++ b/src/openfermion/utils/_sparse_tools.py @@ -598,7 +598,7 @@ def expectation_three_body_db_operator_computational_basis_state( dual_basis_action[4][0] % 2 == orbital3 % 2)): expectation_value += numpy.exp(-1j * ( - k1ad + k2bf + k3ce)) + k1ad + k2bf + k3ce)) # Handle -\delta_{ad} \delta_{be} \delta_{cf} after FT. # The Fourier transform is spin-conserving. From d2630ef6ce3f0332cce9a51c568646a07964a26d Mon Sep 17 00:00:00 2001 From: Ryan Babbush Date: Thu, 28 Sep 2017 19:35:25 -0700 Subject: [PATCH 3/4] addresses the python3 issue with floor division --- src/openfermion/utils/_rdm_equality_constraints.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/openfermion/utils/_rdm_equality_constraints.py b/src/openfermion/utils/_rdm_equality_constraints.py index 2291aa188..847ab32a3 100644 --- a/src/openfermion/utils/_rdm_equality_constraints.py +++ b/src/openfermion/utils/_rdm_equality_constraints.py @@ -76,7 +76,7 @@ def two_body_fermion_constraints(n_orbitals, n_fermions): for ij in range(n_orbitals ** 2): i, j = (ij / n_orbitals), (ij % n_orbitals) for kl in range(ij + 1, n_orbitals ** 2): - k, l = (kl / n_orbitals), (kl % n_orbitals) + k, l = (kl // n_orbitals), (kl % n_orbitals) constraint_operator = FermionOperator( ((i, 1), (j, 1), (l, 0), (k, 0))) constraint_operator -= FermionOperator( @@ -98,9 +98,9 @@ def two_body_fermion_constraints(n_orbitals, n_fermions): # Linear relations between two-particle matrices. for ij in range(n_orbitals ** 2): - i, j = (ij / n_orbitals), (ij % n_orbitals) + i, j = (ij // n_orbitals), (ij % n_orbitals) for kl in range(ij, n_orbitals ** 2): - k, l = (kl / n_orbitals), (kl % n_orbitals) + k, l = (kl // n_orbitals), (kl % n_orbitals) # G-matrix condition. constraint_operator = FermionOperator( From 2edbc67aafcd8da68a0105dc6969ea017ae12105 Mon Sep 17 00:00:00 2001 From: Ryan Babbush Date: Thu, 28 Sep 2017 20:08:10 -0700 Subject: [PATCH 4/4] again, fixing floor divide for python3 --- src/openfermion/utils/_rdm_equality_constraints.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/openfermion/utils/_rdm_equality_constraints.py b/src/openfermion/utils/_rdm_equality_constraints.py index 847ab32a3..de5ad21e2 100644 --- a/src/openfermion/utils/_rdm_equality_constraints.py +++ b/src/openfermion/utils/_rdm_equality_constraints.py @@ -74,7 +74,7 @@ def two_body_fermion_constraints(n_orbitals, n_fermions): # Two-body Hermiticity condition. for ij in range(n_orbitals ** 2): - i, j = (ij / n_orbitals), (ij % n_orbitals) + i, j = (ij // n_orbitals), (ij % n_orbitals) for kl in range(ij + 1, n_orbitals ** 2): k, l = (kl // n_orbitals), (kl % n_orbitals) constraint_operator = FermionOperator(