diff --git a/src/openfermion/circuits/primitives/bogoliubov_transform.py b/src/openfermion/circuits/primitives/bogoliubov_transform.py index be6e6d8b8..ff087308e 100644 --- a/src/openfermion/circuits/primitives/bogoliubov_transform.py +++ b/src/openfermion/circuits/primitives/bogoliubov_transform.py @@ -92,13 +92,14 @@ def bogoliubov_transform( initial_state = _occupied_orbitals(initial_state, n_qubits) initially_occupied_orbitals = cast(Optional[Sequence[int]], initial_state) - # If the transformation matrix is block diagonal with two blocks, - # do each block separately - if _is_spin_block_diagonal(transformation_matrix): - up_block = transformation_matrix[: n_qubits // 2, : n_qubits // 2] + # If the transformation matrix does not mix the two spin sectors, do each + # sector separately. This is both faster (two half-size decompositions) and + # more robust: the general fermionic Gaussian decomposition is numerically + # unreliable on exactly block-structured matrices. See issue #776. + spin_blocks = _spin_blocks(transformation_matrix) + if spin_blocks is not None: + up_block, down_block = spin_blocks up_qubits = qubits[: n_qubits // 2] - - down_block = transformation_matrix[n_qubits // 2 :, n_qubits // 2 :] down_qubits = qubits[n_qubits // 2 :] if initially_occupied_orbitals is None: @@ -110,6 +111,17 @@ def bogoliubov_transform( i - n_qubits // 2 for i in initially_occupied_orbitals if i >= n_qubits // 2 ] + # Spin-down ladder operators carry a Jordan-Wigner Z string over all + # spin-up qubits, so they anticommute with the spin-up parity operator. + # If the spin-up circuit flips spin-up parity (an odd number of + # particle-hole transformations) the factorized circuit would realize + # the spin-down operators with a flipped sign. A layer of Z gates on the + # spin-down qubits restores the documented action U a_p^dagger U^-1 = + # b_p^dagger. When an initial computational basis state is fixed this + # only contributes an irrelevant global phase, so it is skipped there. + if initially_occupied_orbitals is None and not _preserves_parity(up_block): + yield (cirq.Z(q) for q in down_qubits) + yield bogoliubov_transform(up_qubits, up_block, initial_state=up_orbitals) yield bogoliubov_transform(down_qubits, down_block, initial_state=down_orbitals) return @@ -122,6 +134,38 @@ def bogoliubov_transform( yield _gaussian_basis_change(qubits, transformation_matrix, initially_occupied_orbitals) +def _spin_blocks(matrix: numpy.ndarray) -> Optional[Tuple[numpy.ndarray, numpy.ndarray]]: + r"""Split a transformation matrix into independent spin-sector transforms. + + Returns a pair of transformation matrices (spin-up block, spin-down block) + if the transformation does not mix the two spin sectors, assuming modes are + ordered with all spin-up modes before all spin-down modes. Returns None if + the transformation mixes the sectors. + + For an $N \times N$ matrix the blocks are the diagonal blocks. An + $N \times 2N$ matrix $W = (W_1 \; W_2)$ acts on the column vector of all + creation operators followed by all annihilation operators, so it preserves + the spin sectors if and only if both $W_1$ and $W_2$ are block diagonal; + the sector transforms are then $(N/2) \times N$ matrices formed from the + corresponding diagonal blocks of $W_1$ and $W_2$. + """ + n, m = matrix.shape + if n % 2: + return None + k = n // 2 + if m == n: + if _is_spin_block_diagonal(matrix): + return matrix[:k, :k], matrix[k:, k:] + return None + left_block = matrix[:, :n] + right_block = matrix[:, n:] + if _is_spin_block_diagonal(left_block) and _is_spin_block_diagonal(right_block): + up_block = numpy.concatenate([left_block[:k, :k], right_block[:k, :k]], axis=1) + down_block = numpy.concatenate([left_block[k:, k:], right_block[k:, k:]], axis=1) + return up_block, down_block + return None + + def _is_spin_block_diagonal(matrix) -> bool: n = matrix.shape[0] if n % 2: @@ -131,6 +175,30 @@ def _is_spin_block_diagonal(matrix) -> bool: return numpy.isclose(max_upper_right, 0.0) and numpy.isclose(max_lower_left, 0.0) +def _preserves_parity(transformation_matrix: numpy.ndarray) -> bool: + r"""Whether the Bogoliubov circuit for this block preserves fermion parity. + + A particle-number-conserving (square) transformation uses only Givens + rotations and always preserves parity. A general $N \times 2N$ Gaussian + transformation preserves parity if and only if the orthogonal matrix that + implements it on the Majorana operators lies in the identity component of + O(2N), i.e. has determinant $+1$; the other component requires an odd + number of particle-hole transformations, each an X gate that flips the + parity. That determinant equals the determinant of the unitary + Bogoliubov-de Gennes extension of $W = (W_1 \; W_2)$, namely + $\begin{pmatrix} W_1 & W_2 \\ W_2^* & W_1^* \end{pmatrix}$, which is real and + equal to $\pm 1$, so the parity is read off directly without running the + full fermionic Gaussian decomposition. + """ + n, m = transformation_matrix.shape + if m == n: + return True + w1 = transformation_matrix[:, :n] + w2 = transformation_matrix[:, n:] + bogoliubov_de_gennes = numpy.block([[w1, w2], [numpy.conjugate(w2), numpy.conjugate(w1)]]) + return numpy.isclose(numpy.linalg.det(bogoliubov_de_gennes), 1.0) + + def _occupied_orbitals(computational_basis_state: int, n_qubits) -> List[int]: """Indices of ones in the binary expansion of an integer in big endian order. e.g. 010110 -> [1, 3, 4]""" diff --git a/src/openfermion/circuits/primitives/bogoliubov_transform_test.py b/src/openfermion/circuits/primitives/bogoliubov_transform_test.py index 3ed36109f..28349045e 100644 --- a/src/openfermion/circuits/primitives/bogoliubov_transform_test.py +++ b/src/openfermion/circuits/primitives/bogoliubov_transform_test.py @@ -116,6 +116,182 @@ def test_spin_symmetric_bogoliubov_transform( numpy.testing.assert_allclose(quad_ham_sparse.dot(state), energy * state, atol=atol) +def _independent_spin_sector_transform(n_spatial_orbitals, seed_up, seed_down, real): + """Assemble an N x 2N Bogoliubov transformation that does not mix spin. + + Builds two independent non-particle-conserving quadratic Hamiltonians, one + per spin sector, and places their diagonalizing transforms into a combined + matrix. The matrix acts on the column vector of all creation operators + followed by all annihilation operators, so each sector's creation + (annihilation) block lands in the left (right) half, offset by the sector's + position within that half. + """ + quad_ham_up = random_quadratic_hamiltonian( + n_spatial_orbitals, conserves_particle_number=False, real=real, seed=seed_up + ) + quad_ham_down = random_quadratic_hamiltonian( + n_spatial_orbitals, conserves_particle_number=False, real=real, seed=seed_down + ) + up_energies, up_matrix, up_constant = quad_ham_up.diagonalizing_bogoliubov_transform() + down_energies, down_matrix, down_constant = quad_ham_down.diagonalizing_bogoliubov_transform() + + n = n_spatial_orbitals + n_qubits = 2 * n + transformation_matrix = numpy.zeros((n_qubits, 2 * n_qubits), dtype=complex) + transformation_matrix[:n, :n] = up_matrix[:, :n] + transformation_matrix[:n, n_qubits : n_qubits + n] = up_matrix[:, n:] + transformation_matrix[n:, n:n_qubits] = down_matrix[:, :n] + transformation_matrix[n:, n_qubits + n :] = down_matrix[:, n:] + return ( + transformation_matrix, + (quad_ham_up, up_energies, up_constant), + (quad_ham_down, down_energies, down_constant), + ) + + +def test_bogoliubov_transform_spin_block_gaussian_regression(): + """Regression test for issue #776. + + A non-square transformation matrix that does not mix spin sectors used to + be misidentified as a square spin-block-diagonal matrix, leading to bogus + slicing and a ValueError when constructing the circuit. + """ + qubits = LineQubit.range(2) + phase = -9.57167901e-01 - 2.89533435e-01j + phase /= abs(phase) + transformation_matrix = numpy.array([[phase, 0, 0, 0], [0, 1, 0, 0]], dtype=complex) + + circuit = cirq.Circuit( + bogoliubov_transform(qubits, transformation_matrix), cirq.I.on_each(*qubits) + ) + + # The transformation is a phase rotation of mode 0, so it should leave + # computational basis states invariant up to phase. + state = circuit.final_state_vector(initial_state=0, qubit_order=qubits, dtype=numpy.complex128) + expected_state = numpy.zeros(4, dtype=numpy.complex128) + expected_state[0] = 1 + cirq.testing.assert_allclose_up_to_global_phase(state, expected_state, atol=1e-6) + + +@pytest.mark.parametrize( + 'n_spatial_orbitals, real, up_orbitals, down_orbitals', + [(2, True, [0], [0, 1]), (2, False, [1], []), (3, True, [0, 2], [1]), (3, False, [], [0, 2])], +) +def test_bogoliubov_transform_spin_block_gaussian( + n_spatial_orbitals, real, up_orbitals, down_orbitals, atol=5e-5 +): + """A non-particle-conserving transform that does not mix spin sectors must + prepare eigenstates of the combined Hamiltonian with the correct energy.""" + n_qubits = 2 * n_spatial_orbitals + qubits = LineQubit.range(n_qubits) + sim = cirq.Simulator(dtype=numpy.complex128) + + transformation_matrix, up_data, down_data = _independent_spin_sector_transform( + n_spatial_orbitals, 51624, 48201, real + ) + quad_ham_up, up_energies, up_constant = up_data + quad_ham_down, down_energies, down_constant = down_data + + # Combined Hamiltonian with spin-down modes shifted to come after spin-up + ferm_op = openfermion.get_fermion_operator(quad_ham_up) + for term, coefficient in openfermion.get_fermion_operator(quad_ham_down).terms.items(): + shifted_term = tuple((index + n_spatial_orbitals, action) for index, action in term) + ferm_op += openfermion.FermionOperator(shifted_term, coefficient) + combined_ham_sparse = get_sparse_operator(ferm_op, n_qubits=n_qubits) + + energy = ( + sum(up_energies[i] for i in up_orbitals) + + sum(down_energies[i] for i in down_orbitals) + + up_constant + + down_constant + ) + occupied_orbitals = list(up_orbitals) + [i + n_spatial_orbitals for i in down_orbitals] + initial_state = sum(2 ** (n_qubits - 1 - int(i)) for i in occupied_orbitals) + + circuit = cirq.Circuit( + bogoliubov_transform(qubits, transformation_matrix, initial_state=initial_state) + ) + state = sim.simulate( + circuit, initial_state=initial_state, qubit_order=qubits + ).final_state_vector + + numpy.testing.assert_allclose(combined_ham_sparse.dot(state), energy * state, atol=atol) + + +@pytest.mark.parametrize( + 'n_spatial_orbitals, seed_up, seed_down, real', + [ + # The (2, 1000, ...) case has a parity-preserving spin-up sector while + # the (2, 1001, ...) case has a parity-flipping one; the latter is the + # regression guard for the spin-down Jordan-Wigner sign correction. + (2, 1000, 2000, True), + (2, 1001, 2001, True), + (3, 1000, 2000, False), + (3, 1001, 2001, True), + ], +) +def test_bogoliubov_transform_spin_block_operator_algebra( + n_spatial_orbitals, seed_up, seed_down, real, atol=1e-6 +): + """The factorized circuit must implement the documented transformation + exactly: U a_p^dagger U^-1 = b_p^dagger for every mode p, including the + cross-sector Jordan-Wigner sign carried by spin-down operators. Eigenstate + energies cannot detect a wrong sign on b_p^dagger, so this checks the full + operator algebra of the circuit's unitary directly. + """ + n_qubits = 2 * n_spatial_orbitals + qubits = LineQubit.range(n_qubits) + transformation_matrix, _, _ = _independent_spin_sector_transform( + n_spatial_orbitals, seed_up, seed_down, real + ) + + circuit = cirq.Circuit( + bogoliubov_transform(qubits, transformation_matrix), cirq.I.on_each(*qubits) + ) + unitary = circuit.unitary(qubit_order=qubits) + + def ladder(mode, action): + return get_sparse_operator( + openfermion.FermionOperator(((mode, action),)), n_qubits=n_qubits + ).toarray() + + for p in range(n_qubits): + transformed = unitary @ ladder(p, 1) @ unitary.conj().T + expected = numpy.zeros_like(transformed) + for q in range(n_qubits): + expected += transformation_matrix[p, q] * ladder(q, 1) + expected += transformation_matrix[p, n_qubits + q] * ladder(q, 0) + numpy.testing.assert_allclose(transformed, expected, atol=atol) + + +def test_bogoliubov_transform_spin_mixing_gaussian_not_split(atol=5e-5): + """A transform whose annihilation block mixes spin sectors must use the + general procedure even if its creation block is spin-block-diagonal.""" + n_qubits = 4 + qubits = LineQubit.range(n_qubits) + sim = cirq.Simulator(dtype=numpy.complex128) + + # A pairing Hamiltonian whose antisymmetric part couples the two sectors + quad_ham = random_quadratic_hamiltonian( + n_qubits, conserves_particle_number=False, real=True, seed=63003 + ) + quad_ham_sparse = get_sparse_operator(quad_ham) + energies, transformation_matrix, constant = quad_ham.diagonalizing_bogoliubov_transform() + + occupied_orbitals = [0, 2] + energy = sum(energies[i] for i in occupied_orbitals) + constant + initial_state = sum(2 ** (n_qubits - 1 - int(i)) for i in occupied_orbitals) + + circuit = cirq.Circuit( + bogoliubov_transform(qubits, transformation_matrix, initial_state=initial_state) + ) + state = sim.simulate( + circuit, initial_state=initial_state, qubit_order=qubits + ).final_state_vector + + numpy.testing.assert_allclose(quad_ham_sparse.dot(state), energy * state, atol=atol) + + @pytest.mark.parametrize( 'n_qubits, conserves_particle_number', [(4, True), (4, False), (5, True), (5, False)] )