Repository: quantumlib/OpenFermion

Issue:
Circular dependency when importing from utils in ops
I'm working on the file `ops/_quadratic_hamiltonian.py`. In there, I add the line

```from openfermion.utils._slater_determinants import swap_columns, swap_rows```

This causes the tests in the `hamiltonians` folder to fail; i.e., when I run `pytest hamiltonians/`, I get

```
________________ ERROR collecting src/openfermion/hamiltonians/_chemical_series_test.py ________________
ImportError while importing test module '/home/kjs/Projects/OpenFermion/src/openfermion/hamiltonians/_chemical_series_test.py'.
Hint: make sure your test modules/packages have valid Python names.
Traceback:
hamiltonians/__init__.py:13: in <module>
    from ._chemical_series import (make_atomic_ring,
hamiltonians/_chemical_series.py:18: in <module>
    from openfermion.hamiltonians._molecular_data import (MolecularData,
hamiltonians/_molecular_data.py:21: in <module>
    from openfermion.ops import InteractionOperator, InteractionRDM
ops/__init__.py:23: in <module>
    from ._quadratic_hamiltonian import QuadraticHamiltonian
ops/_quadratic_hamiltonian.py:22: in <module>
    from openfermion.utils._slater_determinants import swap_columns, swap_rows
utils/__init__.py:39: in <module>
    from ._trotter_error import error_bound, error_operator
utils/_trotter_error.py:20: in <module>
    from openfermion.hamiltonians import MolecularData
E   ImportError: cannot import name 'MolecularData'
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Interrupted: 1 errors during collection !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
```

It seems that there is a circular dependency
MolecularData -> ops -> utils -> MolecularData
 and any file in the `ops` folder which imports anything from `utils` will cause this error. What is the best way to deal with it?

Discussion:

Patch:
diff --git a/src/openfermion/ops/_quadratic_hamiltonian.py b/src/openfermion/ops/_quadratic_hamiltonian.py
index 15cc467..fdd42be 100644
--- a/src/openfermion/ops/_quadratic_hamiltonian.py
+++ b/src/openfermion/ops/_quadratic_hamiltonian.py
@@ -17,7 +17,7 @@ from __future__ import absolute_import
 import numpy
 
 from openfermion.config import EQ_TOLERANCE
-from openfermion.ops import PolynomialTensor
+from openfermion.ops import FermionOperator, PolynomialTensor
 
 
 class QuadraticHamiltonianError(Exception):
@@ -67,15 +67,17 @@ class QuadraticHamiltonian(PolynomialTensor):
                     hermitian_part -
                     chemical_potential * numpy.eye(n_qubits))
 
-        # Initialize antisymmetric part
         if antisymmetric_part is None:
-            antisymmetric_part = numpy.zeros((n_qubits, n_qubits), complex)
+            super(QuadraticHamiltonian, self).__init__(
+                    {(): constant,
+                     (1, 0): combined_hermitian_part})
+        else:
+            super(QuadraticHamiltonian, self).__init__(
+                    {(): constant,
+                     (1, 0): combined_hermitian_part,
+                     (1, 1): .5 * antisymmetric_part,
+                     (0, 0): -.5 * antisymmetric_part.conj()})
 
-        super(QuadraticHamiltonian, self).__init__(
-                {(): constant,
-                 (1, 0): combined_hermitian_part,
-                 (1, 1): .5 * antisymmetric_part,
-                 (0, 0): -.5 * antisymmetric_part.conj()})
         self.chemical_potential = chemical_potential
 
     def combined_hermitian_part(self):
@@ -92,9 +94,100 @@ class QuadraticHamiltonian(PolynomialTensor):
 
     def antisymmetric_part(self):
         """Return the antisymmetric part."""
-        return 2. * self.n_body_tensors[1, 1].copy()
+        if (1, 1) in self.n_body_tensors:
+            return 2. * self.n_body_tensors[1, 1].copy()
+        else:
+            return numpy.zeros((self.n_qubits, self.n_qubits), complex)
 
     def conserves_particle_number(self):
         """Return whether this Hamiltonian conserves particle number."""
         discrepancy = numpy.max(numpy.abs(self.antisymmetric_part()))
         return discrepancy < EQ_TOLERANCE
+
+    def majorana_form(self):
+        """Return the Majorana represention of the Hamiltonian.
+
+        Any quadratic Hamiltonian can be written in the form
+
+            constant + i / 2 \sum_{j, k} A_{jk} s_j s_k.
+
+        where the s_i are normalized Majorana fermion operators:
+
+            s_j = i / sqrt(2) (a^\dagger_j - a_j)
+            s_{j + n_qubits} = 1 / sqrt(2) (a^\dagger_j + a_j)
+
+        and A is a (2 * n_qubits) x (2 * n_qubits) real antisymmetric matrix.
+        This function returns the matrix A and the constant.
+        """
+        hermitian_part = self.combined_hermitian_part()
+        antisymmetric_part = self.antisymmetric_part()
+
+        # Compute the Majorana matrix using block matrix manipulations
+        majorana_matrix = numpy.zeros((2 * self.n_qubits, 2 * self.n_qubits))
+        # Set upper left block
+        majorana_matrix[:self.n_qubits, :self.n_qubits] = numpy.real(-.5j * (
+                hermitian_part - hermitian_part.conj() +
+                antisymmetric_part - antisymmetric_part.conj()))
+        # Set upper right block
+        majorana_matrix[:self.n_qubits, self.n_qubits:] = numpy.real(.5 * (
+                hermitian_part + hermitian_part.conj() -
+                antisymmetric_part - antisymmetric_part.conj()))
+        # Set lower left block
+        majorana_matrix[self.n_qubits:, :self.n_qubits] = numpy.real(-.5 * (
+                hermitian_part + hermitian_part.conj() +
+                antisymmetric_part + antisymmetric_part.conj()))
+        # Set lower right block
+        majorana_matrix[self.n_qubits:, self.n_qubits:] = numpy.real(-.5j * (
+                hermitian_part - hermitian_part.conj() -
+                antisymmetric_part + antisymmetric_part.conj()))
+
+        # Compute the constant
+        majorana_constant = (.5 * numpy.real(numpy.trace(hermitian_part)) +
+                             self.n_body_tensors[()])
+
+        return majorana_matrix, majorana_constant
+
+
+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.
+
+    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.')
diff --git a/src/openfermion/utils/__init__.py b/src/openfermion/utils/__init__.py
index 9d5bb20..d75656f 100644
--- a/src/openfermion/utils/__init__.py
+++ b/src/openfermion/utils/__init__.py
@@ -20,8 +20,9 @@ from ._operator_utils import (commutator, count_qubits,
                               get_file_path, inverse_fourier_transform,
                               is_identity, load_operator, save_operator)
 
-from ._slater_determinants import (givens_decomposition,
-                                   fermionic_gaussian_decomposition)
+from ._slater_determinants import (fermionic_gaussian_decomposition,
+                                   givens_decomposition,
+                                   ground_state_preparation_circuit)
 
 from ._sparse_tools import (expectation,
                             expectation_computational_basis_state,
diff --git a/src/openfermion/utils/_slater_determinants.py b/src/openfermion/utils/_slater_determinants.py
index 5ebecd6..b3ae971 100644
--- a/src/openfermion/utils/_slater_determinants.py
+++ b/src/openfermion/utils/_slater_determinants.py
@@ -18,6 +18,63 @@ import numpy
 from scipy.linalg import schur
 
 from openfermion.config import EQ_TOLERANCE
+from openfermion.ops import QuadraticHamiltonian
+
+
+def ground_state_preparation_circuit(quadratic_hamiltonian):
+    """Obtain a description of a circuit which prepares the ground state
+    of a quadratic Hamiltonian.
+
+    Args:
+        quadratic_hamiltonian(QuadraticHamiltonian):
+            The Hamiltonian whose ground state is desired.
+
+    Returns:
+        circuit_description(list[tuple]):
+            A list of operations describing the circuit. Each operation
+            is a tuple of objects describing elementary operations that
+            can be performed in parallel. Each elementary operation
+            is either the string 'pht', indicating a particle-hole
+            transformation on the last fermionic mode, or a tuple of
+            the form (i, j, theta, phi), indicating a Givens rotation
+            of qubits i and j by angles theta and phi.
+    """
+    if not isinstance(quadratic_hamiltonian, QuadraticHamiltonian):
+        raise ValueError('Input must be an instance of QuadraticHamiltonian.')
+
+    if quadratic_hamiltonian.conserves_particle_number():
+        # The Hamiltonian conserves particle number, so we don't need
+        # to use the most general procedure.
+        hermitian_matrix = quadratic_hamiltonian.combined_hermitian_part()
+        # Get the unitary rows which represent the ground state
+        # Slater determinant
+        energies, diagonalizing_unitary = numpy.linalg.eigh(
+                hermitian_matrix)
+        num_negative_energies = numpy.count_nonzero(
+                energies < -EQ_TOLERANCE)
+        slater_determinant_matrix = diagonalizing_unitary.T[
+                :num_negative_energies]
+        # Get the circuit description
+        left_unitary, decomposition, diagonal = givens_decomposition(
+                slater_determinant_matrix)
+        circuit_description = decomposition
+    else:
+        # The Hamiltonian does not conserve particle number, so we
+        # need to use the most general procedure.
+        majorana_matrix, majorana_constant = (
+                quadratic_hamiltonian.majorana_form())
+        # Get the unitary rows which represent the ground state
+        # Slater determinant
+        diagonalizing_unitary = diagonalizing_fermionic_unitary(
+                majorana_matrix)
+        slater_determinant_matrix = diagonalizing_unitary[
+                quadratic_hamiltonian.n_qubits:]
+        # Get the circuit description
+        left_unitary, decomposition, diagonal = (
+                fermionic_gaussian_decomposition(
+                    slater_determinant_matrix))
+        circuit_description = decomposition
+    return circuit_description
 
 
 def fermionic_gaussian_decomposition(unitary_rows):
@@ -49,7 +106,7 @@ def fermionic_gaussian_decomposition(unitary_rows):
     something like [('pht', ), (G_1, ), ('pht', G_2), ... ].
     The objects within a tuple are either the string 'pht', which indicates
     a particle-hole transformation on the last fermionic mode, or a tuple
-    of the form (i, j, theta, phi), which indicates a Givens roation
+    of the form (i, j, theta, phi), which indicates a Givens rotation
     of rows i and j by angles theta and phi.
 
     Args:
@@ -148,109 +205,6 @@ def fermionic_gaussian_decomposition(unitary_rows):
     return left_unitary, decomposition, diagonal
 
 
-def diagonalizing_fermionic_unitary(antisymmetric_matrix):
-    """Compute the unitary that diagonalizes a quadratic Hamiltonian.
-
-    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.
-
-    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
-
-
-def antisymmetric_canonical_form(antisymmetric_matrix):
-    """Compute the canonical form of an antisymmetric matrix.
-
-    The input is a real, antisymmetric n x n matrix A, where n is even.
-    Its canonical form is::
-
-        A = R^T C R
-
-    where R is a real, orthogonal matrix and C has the form::
-
-        [  0     D ]
-        [ -D     0 ]
-
-    where D is a diagonal matrix with nonnegative entries.
-
-    Args:
-        antisymmetric_matrix(ndarray): An antisymmetric matrix with even
-            dimension.
-
-    Returns:
-        canonical(ndarray): The canonical form C of antisymmetric_matrix
-        orthogonal(ndarray): The orthogonal transformation R.
-    """
-    m, n = antisymmetric_matrix.shape
-
-    if m != n or n % 2 != 0:
-        raise ValueError('The input matrix must be square with even '
-                         'dimension.')
-
-    # Check that input matrix is antisymmetric
-    matrix_plus_transpose = antisymmetric_matrix + antisymmetric_matrix.T
-    maxval = numpy.max(numpy.abs(matrix_plus_transpose))
-    if maxval > EQ_TOLERANCE:
-        raise ValueError('The input matrix must be antisymmetric.')
-
-    # Compute Schur decomposition
-    canonical, orthogonal = schur(antisymmetric_matrix, output='real')
-
-    # The returned form is block diagonal; we need to permute rows and columns
-    # to put it into the form we want
-    num_blocks = n // 2
-    for i in range(1, num_blocks, 2):
-        swap_rows(canonical, i, num_blocks + i - 1)
-        swap_columns(canonical, i, num_blocks + i - 1)
-        swap_columns(orthogonal, i, num_blocks + i - 1)
-        if num_blocks % 2 != 0:
-            swap_rows(canonical, num_blocks - 1, num_blocks + i)
-            swap_columns(canonical, num_blocks - 1, num_blocks + i)
-            swap_columns(orthogonal, num_blocks - 1, num_blocks + i)
-
-    # Now we permute so that the upper right block is non-negative
-    for i in range(num_blocks):
-        if canonical[i, num_blocks + i] < -EQ_TOLERANCE:
-            swap_rows(canonical, i, num_blocks + i)
-            swap_columns(canonical, i, num_blocks + i)
-            swap_columns(orthogonal, i, num_blocks + i)
-
-    return canonical, orthogonal.T
-
-
 def givens_decomposition(unitary_rows):
     """Decompose a matrix into a sequence of Givens rotations.
 
@@ -379,6 +333,109 @@ def givens_decomposition(unitary_rows):
     return left_unitary, givens_rotations, diagonal
 
 
+def diagonalizing_fermionic_unitary(antisymmetric_matrix):
+    """Compute the unitary that diagonalizes a quadratic Hamiltonian.
+
+    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.
+
+    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
+
+
+def antisymmetric_canonical_form(antisymmetric_matrix):
+    """Compute the canonical form of an antisymmetric matrix.
+
+    The input is a real, antisymmetric n x n matrix A, where n is even.
+    Its canonical form is::
+
+        A = R^T C R
+
+    where R is a real, orthogonal matrix and C has the form::
+
+        [  0     D ]
+        [ -D     0 ]
+
+    where D is a diagonal matrix with nonnegative entries.
+
+    Args:
+        antisymmetric_matrix(ndarray): An antisymmetric matrix with even
+            dimension.
+
+    Returns:
+        canonical(ndarray): The canonical form C of antisymmetric_matrix
+        orthogonal(ndarray): The orthogonal transformation R.
+    """
+    m, n = antisymmetric_matrix.shape
+
+    if m != n or n % 2 != 0:
+        raise ValueError('The input matrix must be square with even '
+                         'dimension.')
+
+    # Check that input matrix is antisymmetric
+    matrix_plus_transpose = antisymmetric_matrix + antisymmetric_matrix.T
+    maxval = numpy.max(numpy.abs(matrix_plus_transpose))
+    if maxval > EQ_TOLERANCE:
+        raise ValueError('The input matrix must be antisymmetric.')
+
+    # Compute Schur decomposition
+    canonical, orthogonal = schur(antisymmetric_matrix, output='real')
+
+    # The returned form is block diagonal; we need to permute rows and columns
+    # to put it into the form we want
+    num_blocks = n // 2
+    for i in range(1, num_blocks, 2):
+        swap_rows(canonical, i, num_blocks + i - 1)
+        swap_columns(canonical, i, num_blocks + i - 1)
+        swap_columns(orthogonal, i, num_blocks + i - 1)
+        if num_blocks % 2 != 0:
+            swap_rows(canonical, num_blocks - 1, num_blocks + i)
+            swap_columns(canonical, num_blocks - 1, num_blocks + i)
+            swap_columns(orthogonal, num_blocks - 1, num_blocks + i)
+
+    # Now we permute so that the upper right block is non-negative
+    for i in range(num_blocks):
+        if canonical[i, num_blocks + i] < -EQ_TOLERANCE:
+            swap_rows(canonical, i, num_blocks + i)
+            swap_columns(canonical, i, num_blocks + i)
+            swap_columns(orthogonal, i, num_blocks + i)
+
+    return canonical, orthogonal.T
+
+
 def givens_matrix_elements(a, b):
     """Compute the matrix elements of the Givens rotation that zeroes out one
     of two row entries.


Test patch:
diff --git a/src/openfermion/ops/_quadratic_hamiltonian_test.py b/src/openfermion/ops/_quadratic_hamiltonian_test.py
index e852031..321f925 100644
--- a/src/openfermion/ops/_quadratic_hamiltonian_test.py
+++ b/src/openfermion/ops/_quadratic_hamiltonian_test.py
@@ -16,20 +16,29 @@ from __future__ import absolute_import
 import numpy
 import unittest
 
-from openfermion.ops import QuadraticHamiltonian
+from openfermion.config import EQ_TOLERANCE
+from openfermion.ops import (FermionOperator,
+                             QuadraticHamiltonian,
+                             normal_ordered)
+from openfermion.ops._quadratic_hamiltonian import majorana_operator
+from openfermion.transforms import get_fermion_operator
 
 
 class QuadraticHamiltoniansTest(unittest.TestCase):
 
     def setUp(self):
         self.n_qubits = 5
-        self.constant = 1.
+        self.constant = 1.7
         self.chemical_potential = 2.
 
         # Obtain random Hermitian and antisymmetric matrices
-        rand_mat = numpy.random.randn(self.n_qubits, self.n_qubits)
+        rand_mat_A = numpy.random.randn(self.n_qubits, self.n_qubits)
+        rand_mat_B = numpy.random.randn(self.n_qubits, self.n_qubits)
+        rand_mat = rand_mat_A + 1.j * rand_mat_B
         self.hermitian_mat = rand_mat + rand_mat.T.conj()
-        rand_mat = numpy.random.randn(self.n_qubits, self.n_qubits)
+        rand_mat_A = numpy.random.randn(self.n_qubits, self.n_qubits)
+        rand_mat_B = numpy.random.randn(self.n_qubits, self.n_qubits)
+        rand_mat = rand_mat_A + 1.j * rand_mat_B
         self.antisymmetric_mat = rand_mat - rand_mat.T
 
         self.combined_hermitian = (
@@ -82,3 +91,43 @@ class QuadraticHamiltoniansTest(unittest.TestCase):
         """Test checking whether Hamiltonian conserves particle number."""
         self.assertTrue(self.quad_ham_pc.conserves_particle_number())
         self.assertFalse(self.quad_ham_npc.conserves_particle_number())
+
+    def test_majorana_form(self):
+        """Test getting the Majorana form."""
+        majorana_matrix, majorana_constant = self.quad_ham_npc.majorana_form()
+        # Convert the Majorana form to a FermionOperator
+        majorana_op = FermionOperator((), majorana_constant)
+        for i in range(2 * self.n_qubits):
+            if i < self.n_qubits:
+                left_op = majorana_operator((i, 1))
+            else:
+                left_op = majorana_operator((i - self.n_qubits, 0))
+            for j in range(2 * self.n_qubits):
+                if j < self.n_qubits:
+                    right_op = majorana_operator((j, 1), majorana_matrix[i, j])
+                else:
+                    right_op = majorana_operator((j - self.n_qubits, 0),
+                                                 majorana_matrix[i, j])
+                majorana_op += .5j * left_op * right_op
+        # Get FermionOperator for original Hamiltonian
+        fermion_operator = normal_ordered(
+                get_fermion_operator(self.quad_ham_npc))
+        self.assertTrue(
+                normal_ordered(majorana_op).isclose(fermion_operator))
+
+
+class MajoranaOperatorTest(unittest.TestCase):
+
+    def test_none_term(self):
+        majorana_op = majorana_operator()
+        self.assertTrue(majorana_operator().isclose(FermionOperator()))
+
+    def test_bad_coefficient(self):
+        with self.assertRaises(ValueError):
+            majorana_op = majorana_operator((1, 1), 'a')
+
+    def test_bad_term(self):
+        with self.assertRaises(ValueError):
+            majorana_op = majorana_operator((2, 2))
+        with self.assertRaises(ValueError):
+            majorana_op = majorana_operator('a')
diff --git a/src/openfermion/utils/_slater_determinants_test.py b/src/openfermion/utils/_slater_determinants_test.py
index d0c838c..da94ad1 100644
--- a/src/openfermion/utils/_slater_determinants_test.py
+++ b/src/openfermion/utils/_slater_determinants_test.py
@@ -18,8 +18,10 @@ import unittest
 from scipy.linalg import qr
 
 from openfermion.config import EQ_TOLERANCE
+from openfermion.ops import QuadraticHamiltonian
 from openfermion.utils import (fermionic_gaussian_decomposition,
-                               givens_decomposition)
+                               givens_decomposition,
+                               ground_state_preparation_circuit)
 from openfermion.utils._slater_determinants import (
         antisymmetric_canonical_form,
         diagonalizing_fermionic_unitary,
@@ -28,6 +30,46 @@ from openfermion.utils._slater_determinants import (
         swap_rows)
 
 
+class GroundStatePreparationCircuitTest(unittest.TestCase):
+
+    def setUp(self):
+        self.n_qubits = 5
+        self.constant = 1.7
+        self.chemical_potential = 2.
+
+        # Obtain random Hermitian and antisymmetric matrices
+        rand_mat_A = numpy.random.randn(self.n_qubits, self.n_qubits)
+        rand_mat_B = numpy.random.randn(self.n_qubits, self.n_qubits)
+        rand_mat = rand_mat_A + 1.j * rand_mat_B
+        self.hermitian_mat = rand_mat + rand_mat.T.conj()
+        rand_mat_A = numpy.random.randn(self.n_qubits, self.n_qubits)
+        rand_mat_B = numpy.random.randn(self.n_qubits, self.n_qubits)
+        rand_mat = rand_mat_A + 1.j * rand_mat_B
+        self.antisymmetric_mat = rand_mat - rand_mat.T
+
+        self.combined_hermitian = (
+                self.hermitian_mat -
+                self.chemical_potential * numpy.eye(self.n_qubits))
+
+        # Initialize a particle-number-conserving Hamiltonian
+        self.quad_ham_pc = QuadraticHamiltonian(
+                self.constant, self.hermitian_mat)
+
+        # Initialize a non-particle-number-conserving Hamiltonian
+        self.quad_ham_npc = QuadraticHamiltonian(
+                self.constant, self.hermitian_mat, self.antisymmetric_mat,
+                self.chemical_potential)
+
+    def test_no_error(self):
+        """Test that the procedure runs without error."""
+        # Test a particle-number-conserving Hamiltonian
+        circuit_description = (
+                ground_state_preparation_circuit(self.quad_ham_pc))
+        # Test a non-particle-number-conserving Hamiltonian
+        circuit_description = (
+                ground_state_preparation_circuit(self.quad_ham_npc))
+
+
 class GivensDecompositionTest(unittest.TestCase):
 
     def test_bad_dimensions(self):