Repository: quantumlib/OpenFermion
Pull request #31 (2017-10-13 01:15:50)

Givens rotations
Zhang and I worked on this together, so I want to make sure he gets credit. This is helpful for compiling circuits to prepare Slater determinants. Since it's a little involved, I made a Jupyter notebook to explain what the code is doing. In short, it takes a matrix of orthonormal rows and decomposes it into a unitary times a diagonal matrix times a sequence of Givens rotations, keeping track of which rotations can be done in parallel. The unitary is introduced because it saves some work in preparing the Slater determinant.

The test cases generate a random starting unitary each time using Scipy's QR decomposition algorithm. I'm not sure if that's appropriate; I can also hard-code some matrices in.

I understand that this is kind of a stand-alone contribution right now, but I think it could fit in. For instance, it can be used in conjunction with a routine that diagonalizes InteractionOperators that have only single-body interactions. Zhang and I are also currently working on code to prepare ground states of general quadratic Hamiltonians (not necessarily particle-number-conserving). I think functionalities for dealing with these types of Hamiltonians would be useful.

kevinsung: Ok I changed some of the variable names in the function (I didn't change the tests but let me know if you want that too). However, I left the single capital letter names for matrices in the function description; let me know if you think that should be changed too.
googlebot: So there's good news and bad news.

:thumbsup: The good news is that everyone that needs to sign a CLA (the pull request submitter and all commit authors) have done so.  Everything is all good there.

:confused: The bad news is that it appears that one or more commits were authored by someone other than the pull request submitter.  We need to confirm that they're okay with their commits being contributed to this project.  Please have them confirm that here in the pull request.

*Note to project maintainer: This is a terminal state, meaning the `cla/google` commit status will not change from this State. It's up to you to confirm consent of the commit author(s) and merge this pull request when appropriate.*

<!-- need_author_consent -->
babbush: overriding CLA who has again lost its mind.

Patch:
diff --git a/examples/givens_rotations.ipynb b/examples/givens_rotations.ipynb
new file mode 100644
index 0000000..919db02
--- /dev/null
+++ b/examples/givens_rotations.ipynb
@@ -0,0 +1,309 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {
+    "collapsed": true
+   },
+   "outputs": [],
+   "source": [
+    "# This code block is for automatic testing purposes, please ignore.\n",
+    "try:\n",
+    "    import openfermion\n",
+    "except:\n",
+    "    import os\n",
+    "    os.chdir('../src/')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Demonstration of Givens rotation decomposition"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {
+    "collapsed": true
+   },
+   "outputs": [],
+   "source": [
+    "import numpy\n",
+    "from openfermion.utils import givens_decomposition\n",
+    "from scipy.linalg import qr\n",
+    "\n",
+    "numpy.set_printoptions(precision=3, linewidth=150, suppress=True)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "A Slater determinant can be represented by an $m \\times n$ matrix with orthonormal rows. Let's obtain such a matrix $Q$ using Scipy's QR decomposition algorithm."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "[[-0.311+0.062j -0.403-0.276j  0.097-0.412j  0.237-0.295j  0.091-0.504j  0.236-0.145j]\n",
+      " [-0.289-0.156j -0.113+0.033j  0.262+0.525j -0.207+0.253j -0.212-0.122j  0.365-0.484j]\n",
+      " [-0.330-0.72j  -0.018-0.099j -0.236-0.182j -0.185+0.005j -0.015+0.025j -0.460-0.164j]]\n"
+     ]
+    }
+   ],
+   "source": [
+    "m, n = (3, 6)\n",
+    "\n",
+    "# Obtain a random matrix of orthonormal rows\n",
+    "x = numpy.random.randn(n, n)\n",
+    "y = numpy.random.randn(n, n)\n",
+    "A = x + 1j*y\n",
+    "Q, R = qr(A)\n",
+    "Q = Q[:m, :]\n",
+    "print(Q)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "There exist unitary matrices $V$ and $U$ such that $VQU^\\dagger$ is an $m \\times n$ diagonal matrix. The matrix $VQ$ represents the same Slater determinant as $Q$ up to an overall phase (equal to $\\det(V)$), but it has zeros in the upper diagonal, which saves some work in the preparation of the Slater determinant. The unitary $U$ can be written in the form\n",
+    "\n",
+    "$$U = G_{N_G} \\cdots G_2 G_1$$\n",
+    "\n",
+    "where the $G_k$ are complex Givens rotations. We can obtain $V$, the $G_k$, and the diagonal entries of $VQU^\\dagger$ using the `givens_decomposition` function."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "[[ 0.240-0.126j -0.439-0.47j  -0.689+0.195j]\n",
+      " [ 0.900+0.065j -0.208+0.072j  0.357-0.101j]\n",
+      " [ 0.335+0.j     0.683+0.268j -0.371+0.461j]]\n",
+      "\n",
+      "[ 0.457+0.889j -0.916-0.4j    0.970-0.242j]\n"
+     ]
+    }
+   ],
+   "source": [
+    "# Get Givens decomposition of U\n",
+    "V, givens_rotations, diagonal = givens_decomposition(Q)\n",
+    "print(V)\n",
+    "print()\n",
+    "print(diagonal)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We didn't print the Givens rotations in the above block because it wouldn't look pretty. The Givens rotations are returned as a list of tuples of tuples. We will now iterate through the tuples of the list, and print the innermost tuple within each tuple as a string."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "('2, 3, 0.671, -2.644',)\n",
+      "('1, 2, 1.421, 0.281', '3, 4, 1.401, -1.617')\n",
+      "('0, 1, 0.683, -2.284', '2, 3, 0.704, 2.724', '4, 5, 1.204, 2.049')\n",
+      "('1, 2, 0.800, -2.721', '3, 4, 1.412, -2.494')\n",
+      "('2, 3, 1.112, -3.003',)\n"
+     ]
+    }
+   ],
+   "source": [
+    "for parallel_set in givens_rotations:\n",
+    "    print(tuple([\"{}, {}, {:.3f}, {:.3f}\".format(i, j, theta, phi)\n",
+    "                 for i, j, theta, phi in parallel_set]))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "There are 5 tuples printed, and within each tuple, there are strings of the form 'i, j, theta, phi'. Each such string represents an innermost tuple, and it is a description of a complex Givens rotation of the coordinates $i$ and $j$ by angles $\\theta$ and $\\varphi$. The $2 \\times 2$ matrix corresponding to this rotation is\n",
+    "\n",
+    "$$\n",
+    "\\begin{pmatrix}\n",
+    "\\cos \\theta & -e^{i \\varphi} \\sin \\theta \\\\\n",
+    "\\sin \\theta & e^{i \\varphi} \\cos \\theta\n",
+    "\\end{pmatrix}.\n",
+    "$$\n",
+    "\n",
+    "The fact that there are 5 tuples means that the circuit depth to prepare the Slater determinant corresponding to $Q$ (up to a phase) has depth 5. All of the rotations within the tuple can be performed in parallel; this is possible because the indices to be rotated are disjoint. For instance, in the third step, we can perform three rotations simultaneously, on coordinates $(0, 1)$, $(2, 3)$, and $(4, 5)$."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Let's check that $VQ$ has zeros in the upper right corner."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "[[ 0.355+0.69j  -0.035+0.087j  0.301-0.385j  0.356-0.154j  0.000+0.j     0.000+0.j   ]\n",
+      " [-0.403-0.176j -0.340-0.323j -0.081-0.496j  0.192-0.297j  0.165-0.427j  0.000-0.j   ]\n",
+      " [ 0.194-0.048j -0.169-0.072j  0.242+0.25j  -0.064-0.069j -0.088-0.325j  0.704-0.432j]]\n"
+     ]
+    }
+   ],
+   "source": [
+    "# Check that VQ has zeros in upper right corner\n",
+    "W = V.dot(Q)\n",
+    "print(W)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Now let's check the Givens decomposition. For each set of Givens rotations that can be performed in parallel, we construct the matrices corresponding to the Givens rotations and multiply them together. Then, we multiply $W = VQ$ repeatedly on the right by these matrices and check that the correct elements are zeroed out."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "metadata": {
+    "scrolled": true
+   },
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Number of rotations to perform in parallel: 1\n",
+      "[[ 0.355+0.69j  -0.035+0.087j  0.384-0.491j  0.000-0.j     0.000+0.j     0.000+0.j   ]\n",
+      " [-0.403-0.176j -0.340-0.323j -0.047-0.607j -0.072-0.032j  0.165-0.427j  0.000-0.j   ]\n",
+      " [ 0.194-0.048j -0.169-0.072j  0.134+0.177j  0.220+0.179j -0.088-0.325j  0.704-0.432j]]\n",
+      "\n",
+      "Number of rotations to perform in parallel: 2\n",
+      "[[ 0.355+0.69j  -0.236+0.585j  0.000-0.j     0.000-0.j     0.000-0.j     0.000+0.j   ]\n",
+      " [-0.403-0.176j  0.160+0.516j -0.368-0.405j -0.425-0.187j  0.000+0.j     0.000-0.j   ]\n",
+      " [ 0.194-0.048j -0.201-0.142j -0.140-0.051j -0.287+0.102j  0.272+0.164j  0.704-0.432j]]\n",
+      "\n",
+      "Number of rotations to perform in parallel: 3\n",
+      "[[ 0.457+0.889j -0.000+0.j     0.000-0.j     0.000-0.j     0.000+0.j     0.000-0.j   ]\n",
+      " [-0.000-0.j    -0.638-0.279j -0.483-0.531j -0.000-0.j     0.000-0.j     0.000+0.j   ]\n",
+      " [ 0.000+0.j     0.308-0.077j -0.303-0.054j  0.141-0.015j  0.759+0.457j  0.000+0.j   ]]\n",
+      "\n",
+      "Number of rotations to perform in parallel: 2\n",
+      "[[ 0.457+0.889j -0.000-0.j    -0.000+0.j     0.000+0.j    -0.000-0.j     0.000-0.j   ]\n",
+      " [-0.000-0.j    -0.916-0.4j    0.000+0.j     0.000-0.j    -0.000-0.j     0.000+0.j   ]\n",
+      " [ 0.000+0.j    -0.000-0.j     0.429-0.107j  0.892-0.095j  0.000-0.j     0.000+0.j   ]]\n",
+      "\n",
+      "Number of rotations to perform in parallel: 1\n",
+      "[[ 0.457+0.889j -0.000-0.j    -0.000+0.j    -0.000+0.j    -0.000-0.j     0.000-0.j   ]\n",
+      " [-0.000-0.j    -0.916-0.4j    0.000-0.j     0.000+0.j    -0.000-0.j     0.000+0.j   ]\n",
+      " [ 0.000+0.j    -0.000-0.j     0.970-0.242j -0.000-0.j     0.000-0.j     0.000+0.j   ]]\n",
+      "\n"
+     ]
+    }
+   ],
+   "source": [
+    "# Check the Givens decomposition\n",
+    "def expanded_givens(G, i, j, n):\n",
+    "    expanded_G = numpy.eye(n, dtype=complex)\n",
+    "    expanded_G[([i], [j]), (i, j)] = G\n",
+    "    return expanded_G\n",
+    "\n",
+    "U = numpy.eye(n, dtype=complex)\n",
+    "for parallel_set in givens_rotations:\n",
+    "    print(\"Number of rotations to perform in parallel: {}\".format(len(parallel_set)))\n",
+    "    combined_givens = numpy.eye(n, dtype=complex)\n",
+    "    for i, j, theta, phi in parallel_set:\n",
+    "        c = numpy.cos(theta)\n",
+    "        s = numpy.sin(theta)\n",
+    "        phase = numpy.exp(1.j * phi)\n",
+    "        G = numpy.array([[c, -phase * s],\n",
+    "                     [s, phase * c]], dtype=complex)\n",
+    "        expanded_G = expanded_givens(G, i, j, n)\n",
+    "        combined_givens = combined_givens.dot(expanded_G)\n",
+    "    W = W.dot(combined_givens.T.conj())\n",
+    "    U = combined_givens.dot(U)\n",
+    "    print(W)\n",
+    "    print()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Finally, let's check that the final matrix, $VQU^\\dagger$, is indeed diagonal, and that its entries match the ones returned by the function."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "metadata": {
+    "scrolled": true
+   },
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "V * Q * U^\\dagger matches the returned diagonal:\n",
+      "True\n"
+     ]
+    }
+   ],
+   "source": [
+    "# Check the diagonal entries\n",
+    "D = numpy.zeros((m, n), dtype=complex)\n",
+    "D[numpy.diag_indices(m)] = diagonal\n",
+    "print(\"V * Q * U^\\dagger matches the returned diagonal:\")\n",
+    "print(numpy.all(numpy.isclose(D, V.dot(Q.dot(U.T.conj())))))"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.6.2"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/src/openfermion/utils/__init__.py b/src/openfermion/utils/__init__.py
index df67604..9c52bb5 100644
--- a/src/openfermion/utils/__init__.py
+++ b/src/openfermion/utils/__init__.py
@@ -14,6 +14,8 @@ from ._chemical_series import (make_atomic_ring,
                                make_atomic_lattice,
                                make_atom)
 
+from ._givens_rotations import givens_decomposition
+
 from ._grid import Grid
 
 from ._hubbard import fermi_hubbard
diff --git a/src/openfermion/utils/_givens_rotations.py b/src/openfermion/utils/_givens_rotations.py
new file mode 100644
index 0000000..e8aa7bf
--- /dev/null
+++ b/src/openfermion/utils/_givens_rotations.py
@@ -0,0 +1,197 @@
+#   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.
+
+"""Givens rotation decomposition of matrices."""
+from __future__ import absolute_import
+
+import numpy
+
+from openfermion.config import EQ_TOLERANCE
+
+
+def givens_decomposition(unitary_rows):
+    """Decompose a matrix into a sequence of Givens rotations.
+
+    The input is an m x n matrix Q with m <= n. The rows of Q are orthonormal.
+    Q can be decomposed as follows:
+
+        V * Q * U^\dagger = D
+
+    where V and U are unitary matrices, and D is an m x n matrix with the
+    first m columns forming a diagonal matrix and the rest of the columns
+    being zero. Furthermore, we can decompose U as
+
+        U = G_k * ... * G_1
+
+    where G_1, ..., G_k are complex Givens rotations, which are invertible
+    n x n matrices. We describe a complex Givens rotation by the column
+    indices (i, j) that it acts on, plus two angles (theta, phi) that
+    characterize the corresponding 2x2 unitary matrix
+
+        [ cos(theta)    -e^{i phi} sin(theta) ]
+        [ sin(theta)     e^{i phi} cos(theta) ]
+
+    Args:
+        unitary_rows: A numpy array or matrix with orthonormal rows,
+            representing the matrix Q.
+    Returns:
+        left_unitary: An m x m numpy array representing the matrix V.
+        givens_rotations: A list of tuples of objects describing Givens
+            rotations. The list looks something like
+            [(G_1, ), (G_2, G_3), ... ]. The Givens rotations within a tuple
+            can be implemented in parallel. The description of a Givens
+            rotation is itself a tuple of the form (i, j, theta, phi), which
+            represents a Givens rotation of rows i and j by angles theta
+            and phi.
+        diagonal: A list of the nonzero entries of D.
+    """
+    current_matrix = numpy.copy(unitary_rows)
+    m, n = current_matrix.shape
+
+    # Check that m <= n
+    if m > n:
+        raise ValueError('The input m x n matrix must have m <= n')
+
+    # Compute left_unitary using Givens rotations
+    left_unitary = numpy.eye(m, dtype=complex)
+    for k in reversed(range(n - m + 1, n)):
+        # Zero out entries in column k
+        for l in range(m - n + k):
+            # Zero out entry in row l
+            givens_rotation = givens_matrix_elements(current_matrix[l, k],
+                                                     current_matrix[l + 1, k])
+            expanded_givens_rotation = expand_two_by_two(givens_rotation,
+                                                         l, l + 1, m)
+            current_matrix = expanded_givens_rotation.dot(current_matrix)
+            left_unitary = expanded_givens_rotation.dot(left_unitary)
+
+    # Compute the decomposition of unitary_rows into Givens rotations
+    givens_rotations = list()
+    if m != n:
+        # There are n - 1 iterations (the circuit depth is n - 1)
+        for k in range(n - 1):
+            # Get the (row, column) indices of elements to zero out in
+            # parallel.
+            # Get the maximum number of simultaneous rotations that
+            # will be performed
+            max_simul_rotations = min(m, n - m)
+            if k < max_simul_rotations - 1:
+                # There are k + 1 elements to zero out
+                start_row = 0
+                end_row = k + 1
+                start_column = n - m - k
+                end_column = start_column + 2 * (k + 1)
+            elif k > n - 1 - max_simul_rotations:
+                # There are n - 1 - k elements to zero out
+                start_row = m - (n - 1 - k)
+                end_row = m
+                start_column = m - (n - 1 - k) + 1
+                end_column = start_column + 2 * (n - 1 - k)
+            else:
+                # There are max_simul_rotations elements to zero out
+                if max_simul_rotations == m:
+                    start_row = 0
+                    end_row = m
+                    start_column = n - m - k
+                    end_column = start_column + 2 * m
+                else:
+                    start_row = k + 1 - max_simul_rotations
+                    end_row = k + 1
+                    start_column = k + 1 - max_simul_rotations + 1
+                    end_column = start_column + 2 * max_simul_rotations
+
+            row_indices = range(start_row, end_row)
+            column_indices = range(start_column, end_column, 2)
+            indices_to_zero_out = zip(row_indices, column_indices)
+
+            parallel_rotations = list()
+            for i, j in indices_to_zero_out:
+                # Compute the Givens rotation to zero out the (i, j) element
+                left_element = current_matrix[i, j - 1].conj()
+                right_element = current_matrix[i, j].conj()
+                givens_rotation = givens_matrix_elements(left_element,
+                                                         right_element)
+                # Need to switch the rows to zero out right_element
+                # rather than left_element
+                givens_rotation = givens_rotation[(1, 0), :]
+
+                # Add the parameters to the list
+                theta = numpy.arccos(numpy.real(givens_rotation[0, 0]))
+                phi = numpy.angle(givens_rotation[1, 1])
+                parallel_rotations.append((j - 1, j, theta, phi))
+
+                # Update the matrix
+                expanded_givens_rotation = expand_two_by_two(givens_rotation,
+                                                             j - 1, j, n)
+                current_matrix = current_matrix.dot(
+                        expanded_givens_rotation.T.conj())
+
+            # Append the current list of parallel rotations to the list
+            givens_rotations.append(tuple(parallel_rotations))
+
+    # Get the diagonal entries
+    diagonal = current_matrix.diagonal()
+
+    return left_unitary, givens_rotations, diagonal
+
+
+def givens_matrix_elements(a, b):
+    """Compute the matrix elements of the Givens rotation that zeroes out one
+    of two row entries.
+
+    Returns a matrix G such that
+
+        G * [a  b]^T= [0  r]^T
+
+    where r is a complex number.
+
+    Args:
+        a: A complex number representing the upper row entry
+        b: A complex number representing the lower row entry
+    Returns:
+        G: A 2 x 2 numpy array representing the matrix G. The numbers in the
+            first column of G are real.
+    """
+    # Handle case that a is zero
+    if abs(a) < EQ_TOLERANCE:
+        c = 1.
+        s = 0.
+        phase = 1.
+    # Handle case that b is zero and a is nonzero
+    elif abs(b) < EQ_TOLERANCE:
+        c = 0.
+        s = 1.
+        sign_a = a / abs(a)
+        phase = sign_a
+    # Handle case that a and b are both nonzero
+    else:
+        denominator = numpy.sqrt(abs(a) ** 2 + abs(b) ** 2)
+
+        c = abs(b) / denominator
+        s = abs(a) / denominator
+        sign_b = b / abs(b)
+        sign_a = a / abs(a)
+        phase = sign_a * sign_b.conjugate()
+
+    # Construct matrix and return
+    givens_rotation = numpy.array([[c, -phase * s],
+                                  [s, phase * c]], dtype=complex)
+    return givens_rotation
+
+
+def expand_two_by_two(M, i, j, n):
+    """Expand the 2 x 2 matrix M to an n x n matrix acting on coordinates
+    i and j.
+    """
+    expanded_M = numpy.eye(n, dtype=complex)
+    expanded_M[([i], [j]), (i, j)] = M
+    return expanded_M


Test patch:
diff --git a/src/openfermion/tests/_example_test.py b/src/openfermion/tests/_example_test.py
index 5010864..1f048de 100644
--- a/src/openfermion/tests/_example_test.py
+++ b/src/openfermion/tests/_example_test.py
@@ -28,6 +28,7 @@ class ExampleTest(unittest.TestCase):
         string_length = len(THIS_DIRECTORY)
         self.directory = THIS_DIRECTORY[:(string_length - 15)] + 'examples/'
         self.demo_name = 'openfermion_demo.ipynb'
+        self.givens_rotations_demo_name = 'givens_rotations.ipynb'
 
     def test_demo(self):
         # Determine if python 2 or 3 is being used.
@@ -60,6 +61,37 @@ class ExampleTest(unittest.TestCase):
             errors = []
         self.assertEqual(errors, [])
 
+    def test_givens_rotations_demo(self):
+        # Determine if python 2 or 3 is being used.
+        major_version, minor_version = sys.version_info[:2]
+        if major_version == 2 or minor_version == 6:
+            version = str(major_version)
+
+            # Run ipython notebook via nbconvert and collect output.
+            with tempfile.NamedTemporaryFile(suffix='.ipynb') as output_file:
+                args = ['jupyter',
+                        'nbconvert',
+                        '--to',
+                        'notebook',
+                        '--execute',
+                        '--ExecutePreprocessor.timeout=600',
+                        '--ExecutePreprocessor.kernel_name=python{}'.format(
+                            version),
+                        '--output',
+                        output_file.name,
+                        self.directory + self.givens_rotations_demo_name]
+                subprocess.check_call(args)
+                output_file.seek(0)
+                nb = nbformat.read(output_file, nbformat.current_nbformat)
+
+            # Parse output and make sure there are no errors.
+            errors = [output for cell in nb.cells if "outputs" in cell for
+                      output in cell["outputs"] if
+                      output.output_type == "error"]
+        else:
+            errors = []
+        self.assertEqual(errors, [])
+
     def test_performance_benchmarks(self):
 
         # Import performance benchmarks and seed random number generator.
diff --git a/src/openfermion/utils/_givens_rotations_test.py b/src/openfermion/utils/_givens_rotations_test.py
new file mode 100644
index 0000000..60fcba8
--- /dev/null
+++ b/src/openfermion/utils/_givens_rotations_test.py
@@ -0,0 +1,405 @@
+#   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 givens_rotations.py."""
+from __future__ import absolute_import
+
+import numpy
+import unittest
+from scipy.linalg import qr
+
+from openfermion.utils import givens_decomposition
+from openfermion.utils._givens_rotations import expand_two_by_two
+
+
+class GivensDecompositionTest(unittest.TestCase):
+
+    def test_bad_dimensions(self):
+        m, n = (3, 2)
+
+        # Obtain a random matrix of orthonormal rows
+        x = numpy.random.randn(m, m)
+        y = numpy.random.randn(m, m)
+        A = x + 1.j*y
+        Q, R = qr(A)
+        Q = Q[:m, :n]
+
+        with self.assertRaises(ValueError):
+            V, givens_rotations, diagonal = givens_decomposition(Q)
+
+    def test_identity(self):
+        n = 3
+        Q = numpy.eye(n, dtype=complex)
+        V, givens_rotations, diagonal = givens_decomposition(Q)
+
+        # V should be the identity
+        I = numpy.eye(n, dtype=complex)
+        for i in range(n):
+            for j in range(n):
+                self.assertAlmostEqual(V[i, j], I[i, j])
+
+        # There should be no Givens rotations
+        self.assertEqual(givens_rotations, list())
+
+        # The diagonal should be ones
+        for d in diagonal:
+            self.assertAlmostEqual(d, 1.)
+
+    def test_antidiagonal(self):
+        m, n = (3, 3)
+        Q = numpy.zeros((m, n), dtype=complex)
+        Q[0, 2] = 1.
+        Q[1, 1] = 1.
+        Q[2, 0] = 1.
+        V, givens_rotations, diagonal = givens_decomposition(Q)
+
+        # There should be no Givens rotations
+        self.assertEqual(givens_rotations, list())
+
+        # VQ should equal the diagonal
+        VQ = V.dot(Q)
+        D = numpy.zeros((m, n), dtype=complex)
+        D[numpy.diag_indices(m)] = diagonal
+        for i in range(n):
+            for j in range(n):
+                self.assertAlmostEqual(VQ[i, j], D[i, j])
+
+    def test_3_by_3(self):
+        m, n = (3, 3)
+        # Obtain a random matrix of orthonormal rows
+        x = numpy.random.randn(n, n)
+        y = numpy.random.randn(n, n)
+        A = x + 1.j*y
+        Q, R = qr(A)
+        Q = Q[:m, :]
+
+        # Get Givens decomposition of Q
+        V, givens_rotations, diagonal = givens_decomposition(Q)
+
+        # There should be no Givens rotations
+        self.assertEqual(givens_rotations, list())
+
+        # Compute V * Q * U^\dagger
+        W = V.dot(Q)
+
+        # Construct the diagonal matrix
+        D = numpy.zeros((m, n), dtype=complex)
+        D[numpy.diag_indices(m)] = diagonal
+
+        # Assert that W and D are the same
+        for i in range(m):
+            for j in range(n):
+                self.assertAlmostEqual(D[i, j], W[i, j])
+
+    def test_3_by_4(self):
+        m, n = (3, 4)
+        # Obtain a random matrix of orthonormal rows
+        x = numpy.random.randn(n, n)
+        y = numpy.random.randn(n, n)
+        A = x + 1.j*y
+        Q, R = qr(A)
+        Q = Q[:m, :]
+
+        # Get Givens decomposition of Q
+        V, givens_rotations, diagonal = givens_decomposition(Q)
+
+        # Compute U
+        U = numpy.eye(n, dtype=complex)
+        for parallel_set in givens_rotations:
+            combined_givens = numpy.eye(n)
+            for i, j, theta, phi in parallel_set:
+                c = numpy.cos(theta)
+                s = numpy.sin(theta)
+                phase = numpy.exp(1.j * phi)
+                G = numpy.array([[c, -phase * s],
+                                [s, phase * c]], dtype=complex)
+                expanded_G = expand_two_by_two(G, i, j, n)
+                combined_givens = combined_givens.dot(expanded_G)
+            U = combined_givens.dot(U)
+
+        # Compute V * Q * U^\dagger
+        W = V.dot(Q.dot(U.T.conj()))
+
+        # Construct the diagonal matrix
+        D = numpy.zeros((m, n), dtype=complex)
+        D[numpy.diag_indices(m)] = diagonal
+
+        # Assert that W and D are the same
+        for i in range(m):
+            for j in range(n):
+                self.assertAlmostEqual(D[i, j], W[i, j])
+
+    def test_3_by_5(self):
+        m, n = (3, 5)
+        # Obtain a random matrix of orthonormal rows
+        x = numpy.random.randn(n, n)
+        y = numpy.random.randn(n, n)
+        A = x + 1.j*y
+        Q, R = qr(A)
+        Q = Q[:m, :]
+
+        # Get Givens decomposition of Q
+        V, givens_rotations, diagonal = givens_decomposition(Q)
+
+        # Compute U
+        U = numpy.eye(n, dtype=complex)
+        for parallel_set in givens_rotations:
+            combined_givens = numpy.eye(n)
+            for i, j, theta, phi in parallel_set:
+                c = numpy.cos(theta)
+                s = numpy.sin(theta)
+                phase = numpy.exp(1.j * phi)
+                G = numpy.array([[c, -phase * s],
+                                [s, phase * c]], dtype=complex)
+                expanded_G = expand_two_by_two(G, i, j, n)
+                combined_givens = combined_givens.dot(expanded_G)
+            U = combined_givens.dot(U)
+
+        # Compute V * Q * U^\dagger
+        W = V.dot(Q.dot(U.T.conj()))
+
+        # Construct the diagonal matrix
+        D = numpy.zeros((m, n), dtype=complex)
+        D[numpy.diag_indices(m)] = diagonal
+
+        # Assert that W and D are the same
+        for i in range(m):
+            for j in range(n):
+                self.assertAlmostEqual(D[i, j], W[i, j])
+
+    def test_3_by_6(self):
+        m, n = (3, 6)
+        # Obtain a random matrix of orthonormal rows
+        x = numpy.random.randn(n, n)
+        y = numpy.random.randn(n, n)
+        A = x + 1.j*y
+        Q, R = qr(A)
+        Q = Q[:m, :]
+
+        # Get Givens decomposition of Q
+        V, givens_rotations, diagonal = givens_decomposition(Q)
+
+        # Compute U
+        U = numpy.eye(n, dtype=complex)
+        for parallel_set in givens_rotations:
+            combined_givens = numpy.eye(n)
+            for i, j, theta, phi in parallel_set:
+                c = numpy.cos(theta)
+                s = numpy.sin(theta)
+                phase = numpy.exp(1.j * phi)
+                G = numpy.array([[c, -phase * s],
+                                [s, phase * c]], dtype=complex)
+                expanded_G = expand_two_by_two(G, i, j, n)
+                combined_givens = combined_givens.dot(expanded_G)
+            U = combined_givens.dot(U)
+
+        # Compute V * Q * U^\dagger
+        W = V.dot(Q.dot(U.T.conj()))
+
+        # Construct the diagonal matrix
+        D = numpy.zeros((m, n), dtype=complex)
+        D[numpy.diag_indices(m)] = diagonal
+
+        # Assert that W and D are the same
+        for i in range(m):
+            for j in range(n):
+                self.assertAlmostEqual(D[i, j], W[i, j])
+
+    def test_3_by_7(self):
+        m, n = (3, 7)
+        # Obtain a random matrix of orthonormal rows
+        x = numpy.random.randn(n, n)
+        y = numpy.random.randn(n, n)
+        A = x + 1.j*y
+        Q, R = qr(A)
+        Q = Q[:m, :]
+
+        # Get Givens decomposition of Q
+        V, givens_rotations, diagonal = givens_decomposition(Q)
+
+        # Compute U
+        U = numpy.eye(n, dtype=complex)
+        for parallel_set in givens_rotations:
+            combined_givens = numpy.eye(n)
+            for i, j, theta, phi in parallel_set:
+                c = numpy.cos(theta)
+                s = numpy.sin(theta)
+                phase = numpy.exp(1.j * phi)
+                G = numpy.array([[c, -phase * s],
+                                [s, phase * c]], dtype=complex)
+                expanded_G = expand_two_by_two(G, i, j, n)
+                combined_givens = combined_givens.dot(expanded_G)
+            U = combined_givens.dot(U)
+
+        # Compute V * Q * U^\dagger
+        W = V.dot(Q.dot(U.T.conj()))
+
+        # Construct the diagonal matrix
+        D = numpy.zeros((m, n), dtype=complex)
+        D[numpy.diag_indices(m)] = diagonal
+
+        # Assert that W and D are the same
+        for i in range(m):
+            for j in range(n):
+                self.assertAlmostEqual(D[i, j], W[i, j])
+
+    def test_3_by_8(self):
+        m, n = (3, 8)
+        # Obtain a random matrix of orthonormal rows
+        x = numpy.random.randn(n, n)
+        y = numpy.random.randn(n, n)
+        A = x + 1.j*y
+        Q, R = qr(A)
+        Q = Q[:m, :]
+
+        # Get Givens decomposition of Q
+        V, givens_rotations, diagonal = givens_decomposition(Q)
+
+        # Compute U
+        U = numpy.eye(n, dtype=complex)
+        for parallel_set in givens_rotations:
+            combined_givens = numpy.eye(n)
+            for i, j, theta, phi in parallel_set:
+                c = numpy.cos(theta)
+                s = numpy.sin(theta)
+                phase = numpy.exp(1.j * phi)
+                G = numpy.array([[c, -phase * s],
+                                [s, phase * c]], dtype=complex)
+                expanded_G = expand_two_by_two(G, i, j, n)
+                combined_givens = combined_givens.dot(expanded_G)
+            U = combined_givens.dot(U)
+
+        # Compute V * Q * U^\dagger
+        W = V.dot(Q.dot(U.T.conj()))
+
+        # Construct the diagonal matrix
+        D = numpy.zeros((m, n), dtype=complex)
+        D[numpy.diag_indices(m)] = diagonal
+
+        # Assert that W and D are the same
+        for i in range(m):
+            for j in range(n):
+                self.assertAlmostEqual(D[i, j], W[i, j])
+
+    def test_3_by_9(self):
+        m, n = (3, 9)
+        # Obtain a random matrix of orthonormal rows
+        x = numpy.random.randn(n, n)
+        y = numpy.random.randn(n, n)
+        A = x + 1.j*y
+        Q, R = qr(A)
+        Q = Q[:m, :]
+
+        # Get Givens decomposition of Q
+        V, givens_rotations, diagonal = givens_decomposition(Q)
+
+        # Compute U
+        U = numpy.eye(n, dtype=complex)
+        for parallel_set in givens_rotations:
+            combined_givens = numpy.eye(n)
+            for i, j, theta, phi in parallel_set:
+                c = numpy.cos(theta)
+                s = numpy.sin(theta)
+                phase = numpy.exp(1.j * phi)
+                G = numpy.array([[c, -phase * s],
+                                [s, phase * c]], dtype=complex)
+                expanded_G = expand_two_by_two(G, i, j, n)
+                combined_givens = combined_givens.dot(expanded_G)
+            U = combined_givens.dot(U)
+
+        # Compute V * Q * U^\dagger
+        W = V.dot(Q.dot(U.T.conj()))
+
+        # Construct the diagonal matrix
+        D = numpy.zeros((m, n), dtype=complex)
+        D[numpy.diag_indices(m)] = diagonal
+
+        # Assert that W and D are the same
+        for i in range(m):
+            for j in range(n):
+                self.assertAlmostEqual(D[i, j], W[i, j])
+
+    def test_4_by_5(self):
+        m, n = (4, 5)
+        # Obtain a random matrix of orthonormal rows
+        x = numpy.random.randn(n, n)
+        y = numpy.random.randn(n, n)
+        A = x + 1.j*y
+        Q, R = qr(A)
+        Q = Q[:m, :]
+
+        # Get Givens decomposition of Q
+        V, givens_rotations, diagonal = givens_decomposition(Q)
+
+        # Compute U
+        U = numpy.eye(n, dtype=complex)
+        for parallel_set in givens_rotations:
+            combined_givens = numpy.eye(n)
+            for i, j, theta, phi in parallel_set:
+                c = numpy.cos(theta)
+                s = numpy.sin(theta)
+                phase = numpy.exp(1.j * phi)
+                G = numpy.array([[c, -phase * s],
+                                [s, phase * c]], dtype=complex)
+                expanded_G = expand_two_by_two(G, i, j, n)
+                combined_givens = combined_givens.dot(expanded_G)
+            U = combined_givens.dot(U)
+
+        # Compute V * Q * U^\dagger
+        W = V.dot(Q.dot(U.T.conj()))
+
+        # Construct the diagonal matrix
+        D = numpy.zeros((m, n), dtype=complex)
+        D[numpy.diag_indices(m)] = diagonal
+
+        # Assert that W and D are the same
+        for i in range(m):
+            for j in range(n):
+                self.assertAlmostEqual(D[i, j], W[i, j])
+
+    def test_4_by_9(self):
+        m, n = (4, 9)
+        # Obtain a random matrix of orthonormal rows
+        x = numpy.random.randn(n, n)
+        y = numpy.random.randn(n, n)
+        A = x + 1.j*y
+        Q, R = qr(A)
+        Q = Q[:m, :]
+
+        # Get Givens decomposition of Q
+        V, givens_rotations, diagonal = givens_decomposition(Q)
+
+        # Compute U
+        U = numpy.eye(n, dtype=complex)
+        for parallel_set in givens_rotations:
+            combined_givens = numpy.eye(n)
+            for i, j, theta, phi in parallel_set:
+                c = numpy.cos(theta)
+                s = numpy.sin(theta)
+                phase = numpy.exp(1.j * phi)
+                G = numpy.array([[c, -phase * s],
+                                [s, phase * c]], dtype=complex)
+                expanded_G = expand_two_by_two(G, i, j, n)
+                combined_givens = combined_givens.dot(expanded_G)
+            U = combined_givens.dot(U)
+
+        # Compute V * Q * U^\dagger
+        W = V.dot(Q.dot(U.T.conj()))
+
+        # Construct the diagonal matrix
+        D = numpy.zeros((m, n), dtype=complex)
+        D[numpy.diag_indices(m)] = diagonal
+
+        # Assert that W and D are the same
+        for i in range(m):
+            for j in range(n):
+                self.assertAlmostEqual(D[i, j], W[i, j])