Permalink
Browse files

Minor restructuring (#120)

* reorganized diagonalizing Bogoliubov transformation

* pep8 lint

* added my middle initial to NOTICE and README

* move majorana_operator to _fermion_operator.py and expose it

* added option to orbital_energies to enforce nonnegativity

* don't add empty operations to state preparation circuits

* fixed description of get_sparse_operator

* expanded description

* added ground energy method to QuadraticHamiltonian

* added testing method for generating random unitaries

* renamed test function

* fixed Sphinx warning and removed pep8 lint

* moved test function, removed pep8 lint

* added tests/_testing_utils.py

* moved _testing_utils.py to utils

* import absolute_import

* added utils/_testing_utils.py
  • Loading branch information...
kevinsung authored and babbush committed Dec 11, 2017
1 parent 24c160a commit 62c7ccc8bb3867325926d254c86216d99d745078
View
2 NOTICE
@@ -5,7 +5,7 @@ Ryan Babbush (Google)
Jarrod McClean (Google)
Ian Kivlichan (Harvard)
Damian Steiger (ETH Zurich)
Kevin Sung (University of Michigan)
Kevin J. Sung (University of Michigan)
Dave Bacon (Google)
Yudong Cao (Harvard)
View
@@ -117,7 +117,7 @@ Authors
`Jarrod McClean <http://jarrodmcclean.com>`__ (Google),
`Ian Kivlichan <http://aspuru.chem.harvard.edu/ian-kivlichan/>`__ (Harvard),
`Damian Steiger <https://github.com/damiansteiger>`__ (ETH Zurich),
`Kevin Sung <https://github.com/kevinsung>`__ (University of Michigan),
`Kevin J. Sung <https://github.com/kevinsung>`__ (University of Michigan),
`Dave Bacon <https://github.com/dabacon>`__ (Google),
`Yudong Cao <https://github.com/yudongcao>`__ (Harvard),
`Chengyu Dai <https://github.com/jdaaph>`__ (University of Michigan),
@@ -142,7 +142,7 @@ How to cite
===========
When using OpenFermion for research projects, please cite:
Jarrod R. McClean, Ian D. Kivlichan, Damian S. Steiger, Kevin Sung,
Jarrod R. McClean, Ian D. Kivlichan, Damian S. Steiger, Kevin J. Sung,
Yudong Cao, Chengyu Dai, E. Schuyler Fried, Craig Gidney, Thomas Häner,
Vojtĕch Havlíček, Cupjin Huang, Zhang Jiang, Matthew Neeley, Jhonathan Romero,
Nicholas Rubin, Nicolas P. D. Sawaya, Kanav Setia, Sukin Sim, Wei Sun,
@@ -19,58 +19,7 @@
normal_ordered)
from openfermion.transforms import get_fermion_operator, jordan_wigner
from openfermion.utils import jordan_wigner_sparse
def artificial_molecular_operator(n_qubits):
"""Make an artificial random InteractionOperator for testing purposes."""
# Initialize.
constant = numpy.random.randn()
one_body_coefficients = numpy.zeros((n_qubits, n_qubits), float)
two_body_coefficients = numpy.zeros((n_qubits, n_qubits,
n_qubits, n_qubits), float)
# Randomly generate the one-body and two-body integrals.
for p in range(n_qubits):
for q in range(n_qubits):
# One-body terms.
if (p <= p) and (p % 2 == q % 2):
one_body_coefficients[p, q] = numpy.random.randn()
one_body_coefficients[q, p] = one_body_coefficients[p, q]
# Keep looping.
for r in range(n_qubits):
for s in range(n_qubits):
# Skip zero terms.
if (p == q) or (r == s):
continue
# Identify and skip one of the complex conjugates.
if [p, q, r, s] != [s, r, q, p]:
unique_indices = len(set([p, q, r, s]))
# srqp srpq sprq spqr sqpr sqrp
# rsqp rspq rpsq rpqs rqps rqsp.
if unique_indices == 4:
if min(r, s) <= min(p, q):
continue
# qqpp.
elif unique_indices == 2:
if q < p:
continue
# Add the two-body coefficients.
two_body_coefficients[p, q, r, s] = numpy.random.randn()
two_body_coefficients[s, r, q, p] = two_body_coefficients[
p, q, r, s]
# Build the molecular operator and return.
molecular_operator = InteractionOperator(
constant, one_body_coefficients, two_body_coefficients)
return molecular_operator
from openfermion.utils._testing_utils import random_interaction_operator
def benchmark_molecular_operator_jordan_wigner(n_qubits):
@@ -85,7 +34,7 @@ def benchmark_molecular_operator_jordan_wigner(n_qubits):
runtime: The number of seconds required to make the conversion.
"""
# Get an instance of InteractionOperator.
molecular_operator = artificial_molecular_operator(n_qubits)
molecular_operator = random_interaction_operator(n_qubits)
# Convert to a qubit operator.
start = time.time()
@@ -169,7 +118,7 @@ def benchmark_jordan_wigner_sparse(n_qubits):
runtime: The time in seconds that the benchmark took.
"""
# Initialize a random FermionOperator.
molecular_operator = artificial_molecular_operator(n_qubits)
molecular_operator = random_interaction_operator(n_qubits)
fermion_operator = get_fermion_operator(molecular_operator)
# Map to SparseOperator class.
@@ -12,6 +12,7 @@
from ._fermion_operator import (FermionOperator,
hermitian_conjugated,
majorana_operator,
normal_ordered,
number_operator)
from ._qubit_operator import QubitOperator
@@ -31,6 +31,54 @@ def hermitian_conjugated(fermion_operator):
return conjugate_operator
def majorana_operator(term=None, coefficient=1.):
"""Initialize a Majorana operator.
Args:
term(tuple): The first element of the tuple indicates the mode
on which the Majorana operator acts, starting from zero.
The second element of the tuple is an integer, either 1 or 0,
indicating which type of Majorana operator it is:
Type 1: :math:`\\frac{1}{\sqrt{2}} (a^\dagger_p + a_p)`
Type 0: :math:`\\frac{i}{\sqrt{2}} (a^\dagger_p - a_p)`
where the :math:`a^\dagger_p` and :math:`a_p` are the usual
fermionic ladder operators.
Default will result in the zero operator.
coefficient(complex or float, optional): The coefficient of the term.
Default value is 1.0.
Returns:
FermionOperator
"""
if not isinstance(coefficient, (int, float, complex)):
raise ValueError('Coefficient must be scalar.')
if term is None:
# Return zero operator
return FermionOperator()
elif isinstance(term, tuple):
mode, operator_type = term
if operator_type == 1:
majorana_op = FermionOperator(
((mode, 1),), coefficient / numpy.sqrt(2.))
majorana_op += FermionOperator(
((mode, 0),), coefficient / numpy.sqrt(2.))
elif operator_type == 0:
majorana_op = FermionOperator(
((mode, 1),), 1.j * coefficient / numpy.sqrt(2.))
majorana_op -= FermionOperator(
((mode, 0),), 1.j * coefficient / numpy.sqrt(2.))
else:
raise ValueError('Operator specified incorrectly.')
return majorana_op
else:
# Invalid input.
raise ValueError('Operator specified incorrectly.')
def number_operator(n_orbitals, orbital=None, coefficient=1.):
"""Return a number operator.
@@ -406,7 +406,7 @@ def test_rotate_basis_max_order(self):
for order in [27, 28]:
with self.assertRaises(ValueError):
tensor, want_tensor = self.do_rotate_basis_high_order(order)
def do_rotate_basis_high_order(self, order):
key = (1,) * (order // 2) + (0,) * ((order + 1) // 2)
shape = (1,) * order
@@ -417,7 +417,8 @@ def do_rotate_basis_high_order(self, order):
{key: numpy.zeros(shape) + num})
# If order is odd, there are one more 0 than 1 in key
if order % 2 == 1: num *= rotation
if order % 2 == 1:
num *= rotation
want_polynomial_tensor = PolynomialTensor(
{key: numpy.zeros(shape) + num})
@@ -125,28 +125,34 @@ def add_chemical_potential(self, chemical_potential):
numpy.eye(self.n_qubits))
self.chemical_potential += chemical_potential
def orbital_energies(self):
def orbital_energies(self, non_negative=False):
"""Return the orbital energies.
Any quadratic Hamiltonian is unitarily equivalent to a Hamiltonian
of the form
.. math::
\sum_{j} \\varepsilon_j a^\dagger_j a_j + \\text{constant}.
\sum_{j} \\varepsilon_j b^\dagger_j b_j + \\text{constant}.
We call the :math:`\\varepsilon_j` the orbital energies.
The eigenvalues of the Hamiltonian are sums of subsets of the
orbital energies (up to the additive constant).
Args:
non_negative(bool): If True, always return a list of orbital
energies that are non-negative. This option is ignored if
the Hamiltonian does not conserve particle number, in which
case the returned orbital energies are always non-negative.
Returns
-------
orbital_energies(ndarray)
A one-dimensional array containing the :math:`\\varepsilon_j`
constant(float)
The constant
"""
if self.conserves_particle_number:
if self.conserves_particle_number and not non_negative:
hermitian_matrix = self.combined_hermitian_part
orbital_energies, diagonalizing_unitary = numpy.linalg.eigh(
hermitian_matrix)
@@ -162,6 +168,11 @@ def orbital_energies(self):
return orbital_energies, constant
def ground_energy(self):
"""Return the ground energy."""
orbital_energies, constant = self.orbital_energies(non_negative=True)
return constant
def majorana_form(self):
"""Return the Majorana represention of the Hamiltonian.
@@ -211,93 +222,71 @@ def majorana_form(self):
return majorana_matrix, majorana_constant
def diagonalizing_bogoliubov_transform(self):
"""Compute the unitary that diagonalizes a quadratic Hamiltonian.
def majorana_operator(term=None, coefficient=1.):
"""Initialize a Majorana operator.
Args:
term(tuple): The first element of the tuple indicates the mode
on which the Majorana operator acts, starting from zero.
The second element of the tuple is an integer, either 1 or 0,
indicating which type of Majorana operator it is:
type 1: 1 / sqrt(2) (a^\dagger_j + a_j)
type 0: i / sqrt(2) (a^\dagger_j - a_j)
where the a^\dagger_j and a_j are the usual fermionic ladder
operators.
Default will result in the zero operator.
coefficient(complex or float, optional): The coefficient of the term.
Default value is 1.0.
Any quadratic Hamiltonian can be rewritten in the form
Returns:
FermionOperator
"""
if not isinstance(coefficient, (int, float, complex)):
raise ValueError('Coefficient must be scalar.')
if term is None:
# Return zero operator
return FermionOperator()
elif isinstance(term, tuple):
mode, operator_type = term
if operator_type == 1:
majorana_op = FermionOperator(
((mode, 1),), coefficient / numpy.sqrt(2.))
majorana_op += FermionOperator(
((mode, 0),), coefficient / numpy.sqrt(2.))
elif operator_type == 0:
majorana_op = FermionOperator(
((mode, 1),), 1.j * coefficient / numpy.sqrt(2.))
majorana_op -= FermionOperator(
((mode, 0),), 1.j * coefficient / numpy.sqrt(2.))
else:
raise ValueError('Operator specified incorrectly.')
return majorana_op
# Invalid input.
else:
raise ValueError('Operator specified incorrectly.')
.. math::
\sum_{j} \\varepsilon_j b^\dagger_j b_j + \\text{constant},
def diagonalizing_fermionic_unitary(antisymmetric_matrix):
"""Compute the unitary that diagonalizes a quadratic Hamiltonian.
where the :math:`b_j` are a new set fermionic operators
that satisfy the canonical anticommutation relations.
The new fermionic operators are linear combinations of the
original ones:
The input matrix represents a quadratic Hamiltonian in the Majorana basis.
The output matrix is a unitary that represents a transformation (mixing)
of the fermionic ladder operators. We use the convention that the
creation operators are listed before the annihilation operators.
The returned unitary has additional structure which ensures
that the transformed ladder operators also satisfy the fermionic
anticommutation relations.
.. math::
Args:
antisymmetric_matrix(ndarray): A (2 * n_qubits) x (2 * n_qubits)
antisymmetric matrix representing a quadratic Hamiltonian in the
Majorana basis.
Returns:
diagonalizing_unitary(ndarray): A (2 * n_qubits) x (2 * n_qubits)
unitary matrix representing a transformation of the fermionic
ladder operators.
"""
m, n = antisymmetric_matrix.shape
n_qubits = n // 2
# Get the orthogonal transformation that puts antisymmetric_matrix
# into canonical form
canonical, orthogonal = antisymmetric_canonical_form(antisymmetric_matrix)
# Create the matrix that converts between fermionic ladder and
# Majorana bases
normalized_identity = numpy.eye(n_qubits, dtype=complex) / numpy.sqrt(2.)
majorana_basis_change = numpy.eye(
2 * n_qubits, dtype=complex) / numpy.sqrt(2.)
majorana_basis_change[n_qubits:, n_qubits:] *= -1.j
majorana_basis_change[:n_qubits, n_qubits:] = normalized_identity
majorana_basis_change[n_qubits:, :n_qubits] = 1.j * normalized_identity
# Compute the unitary and return
diagonalizing_unitary = majorana_basis_change.T.conj().dot(
orthogonal.dot(majorana_basis_change))
return diagonalizing_unitary
\\begin{pmatrix}
b^\dagger_1 \\\\
\\vdots \\\\
b^\dagger_N \\\\
b_1 \\\\
\\vdots \\\\
b_N
\\end{pmatrix}
= W
\\begin{pmatrix}
a^\dagger_1 \\\\
\\vdots \\\\
a^\dagger_N \\\\
a_1 \\\\
\\vdots \\\\
a_N
\\end{pmatrix},
where :math:`W` is a :math:`2N \\times 2N` unitary matrix.
This method returns the matrix :math:`W`.
Returns:
diagonalizing_unitary (ndarray):
A (2 * `n_qubits`) x (2 * `n_qubits`) matrix representing
the transformation :math:`W` of the fermionic ladder operators.
"""
majorana_matrix, majorana_constant = self.majorana_form()
# Get the orthogonal transformation that puts majorana_matrix
# into canonical form
canonical, orthogonal = antisymmetric_canonical_form(majorana_matrix)
# Create the matrix that converts between fermionic ladder and
# Majorana bases
normalized_identity = (numpy.eye(self.n_qubits, dtype=complex) /
numpy.sqrt(2.))
majorana_basis_change = numpy.eye(
2 * self.n_qubits, dtype=complex) / numpy.sqrt(2.)
majorana_basis_change[self.n_qubits:, self.n_qubits:] *= -1.j
majorana_basis_change[:self.n_qubits,
self.n_qubits:] = normalized_identity
majorana_basis_change[self.n_qubits:,
:self.n_qubits] = 1.j * normalized_identity
# Compute the unitary and return
diagonalizing_unitary = majorana_basis_change.T.conj().dot(
orthogonal.dot(majorana_basis_change))
return diagonalizing_unitary
def antisymmetric_canonical_form(antisymmetric_matrix):
Oops, something went wrong.

0 comments on commit 62c7ccc

Please sign in to comment.