# 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.
"""This module provides functions to interface with scipy.sparse."""
import itertools
from functools import reduce
import numpy.linalg
import numpy
import scipy
import scipy.sparse
import scipy.sparse.linalg
from openfermion.config import EQ_TOLERANCE
from openfermion.ops import (FermionOperator, QuadraticHamiltonian,
QubitOperator, BosonOperator,
QuadOperator, up_index, down_index)
from openfermion.utils import (count_qubits, gaussian_state_preparation_circuit,
is_hermitian,
slater_determinant_preparation_circuit)
# Make global definitions.
identity_csc = scipy.sparse.identity(2, format='csc', dtype=complex)
pauli_x_csc = scipy.sparse.csc_matrix([[0., 1.], [1., 0.]], dtype=complex)
pauli_y_csc = scipy.sparse.csc_matrix([[0., -1.j], [1.j, 0.]], dtype=complex)
pauli_z_csc = scipy.sparse.csc_matrix([[1., 0.], [0., -1.]], dtype=complex)
q_raise_csc = (pauli_x_csc - 1.j * pauli_y_csc) / 2.
q_lower_csc = (pauli_x_csc + 1.j * pauli_y_csc) / 2.
pauli_matrix_map = {'I': identity_csc, 'X': pauli_x_csc,
'Y': pauli_y_csc, 'Z': pauli_z_csc}
def wrapped_kronecker(operator_1, operator_2):
"""Return the Kronecker product of two sparse.csc_matrix operators."""
return scipy.sparse.kron(operator_1, operator_2, 'csc')
def kronecker_operators(*args):
"""Return the Kronecker product of multiple sparse.csc_matrix operators."""
return reduce(wrapped_kronecker, *args)
def jordan_wigner_ladder_sparse(n_qubits, tensor_factor, ladder_type):
r"""Make a matrix representation of a fermion ladder operator.
Operators are mapped as follows:
a_j^\dagger -> Z_0 .. Z_{j-1} (X_j - iY_j) / 2
a_j -> Z_0 .. Z_{j-1} (X_j + iY_j) / 2
Args:
index: This is a nonzero integer. The integer indicates the tensor
factor and the sign indicates raising or lowering.
n_qubits(int): Number qubits in the system Hilbert space.
Returns:
The corresponding Scipy sparse matrix.
"""
parities = tensor_factor * [pauli_z_csc]
identities = [scipy.sparse.identity(
2 ** (n_qubits - tensor_factor - 1), dtype=complex, format='csc')]
if ladder_type:
operator = kronecker_operators(parities + [q_raise_csc] + identities)
else:
operator = kronecker_operators(parities + [q_lower_csc] + identities)
return operator
def jordan_wigner_sparse(fermion_operator, n_qubits=None):
r"""Initialize a Scipy sparse matrix from a FermionOperator.
Operators are mapped as follows:
a_j^\dagger -> Z_0 .. Z_{j-1} (X_j - iY_j) / 2
a_j -> Z_0 .. Z_{j-1} (X_j + iY_j) / 2
Args:
fermion_operator(FermionOperator): instance of the FermionOperator
class.
n_qubits(int): Number of qubits.
Returns:
The corresponding Scipy sparse matrix.
"""
if n_qubits is None:
n_qubits = count_qubits(fermion_operator)
# Create a list of raising and lowering operators for each orbital.
jw_operators = []
for tensor_factor in range(n_qubits):
jw_operators += [(jordan_wigner_ladder_sparse(n_qubits,
tensor_factor,
0),
jordan_wigner_ladder_sparse(n_qubits,
tensor_factor,
1))]
# Construct the Scipy sparse matrix.
n_hilbert = 2 ** n_qubits
values_list = [[]]
row_list = [[]]
column_list = [[]]
for term in fermion_operator.terms:
coefficient = fermion_operator.terms[term]
sparse_matrix = coefficient * scipy.sparse.identity(
2 ** n_qubits, dtype=complex, format='csc')
for ladder_operator in term:
sparse_matrix = sparse_matrix * jw_operators[
ladder_operator[0]][ladder_operator[1]]
if coefficient:
# Extract triplets from sparse_term.
sparse_matrix = sparse_matrix.tocoo(copy=False)
values_list.append(sparse_matrix.data)
(row, column) = sparse_matrix.nonzero()
row_list.append(row)
column_list.append(column)
values_list = numpy.concatenate(values_list)
row_list = numpy.concatenate(row_list)
column_list = numpy.concatenate(column_list)
sparse_operator = scipy.sparse.coo_matrix((
values_list, (row_list, column_list)),
shape=(n_hilbert, n_hilbert)).tocsc(copy=False)
sparse_operator.eliminate_zeros()
return sparse_operator
def qubit_operator_sparse(qubit_operator, n_qubits=None):
"""Initialize a Scipy sparse matrix from a QubitOperator.
Args:
qubit_operator(QubitOperator): instance of the QubitOperator class.
n_qubits (int): Number of qubits.
Returns:
The corresponding Scipy sparse matrix.
"""
if n_qubits is None:
n_qubits = count_qubits(qubit_operator)
if n_qubits < count_qubits(qubit_operator):
raise ValueError('Invalid number of qubits specified.')
# Construct the Scipy sparse matrix.
n_hilbert = 2 ** n_qubits
values_list = [[]]
row_list = [[]]
column_list = [[]]
# Loop through the terms.
for qubit_term in qubit_operator.terms:
tensor_factor = 0
coefficient = qubit_operator.terms[qubit_term]
sparse_operators = [coefficient]
for pauli_operator in qubit_term:
# Grow space for missing identity operators.
if pauli_operator[0] > tensor_factor:
identity_qubits = pauli_operator[0] - tensor_factor
identity = scipy.sparse.identity(
2 ** identity_qubits, dtype=complex, format='csc')
sparse_operators += [identity]
# Add actual operator to the list.
sparse_operators += [pauli_matrix_map[pauli_operator[1]]]
tensor_factor = pauli_operator[0] + 1
# Grow space at end of string unless operator acted on final qubit.
if tensor_factor < n_qubits or not qubit_term:
identity_qubits = n_qubits - tensor_factor
identity = scipy.sparse.identity(
2 ** identity_qubits, dtype=complex, format='csc')
sparse_operators += [identity]
# Extract triplets from sparse_term.
sparse_matrix = kronecker_operators(sparse_operators)
values_list.append(sparse_matrix.tocoo(copy=False).data)
(column, row) = sparse_matrix.nonzero()
column_list.append(column)
row_list.append(row)
# Create sparse operator.
values_list = numpy.concatenate(values_list)
row_list = numpy.concatenate(row_list)
column_list = numpy.concatenate(column_list)
sparse_operator = scipy.sparse.coo_matrix((
values_list, (row_list, column_list)),
shape=(n_hilbert, n_hilbert)).tocsc(copy=False)
sparse_operator.eliminate_zeros()
return sparse_operator
def get_linear_qubit_operator_diagonal(qubit_operator, n_qubits=None):
""" Return a linear operator's diagonal elements.
The main motivation is to use it for Davidson's algorithm, to find out the
lowest n eigenvalues and associated eigenvectors.
Qubit terms with X or Y operators will contribute nothing to the diagonal
elements, while I or Z will contribute a factor of 1 or -1 together with
the coefficient.
Args:
qubit_operator(QubitOperator): A qubit operator.
Returns:
linear_operator_diagonal(numpy.ndarray): The diagonal elements for
LinearQubitOperator(qubit_operator).
"""
if n_qubits is None:
n_qubits = count_qubits(qubit_operator)
if n_qubits < count_qubits(qubit_operator):
raise ValueError('Invalid number of qubits specified.')
n_hilbert = 2 ** n_qubits
zeros_diagonal = numpy.zeros(n_hilbert)
ones_diagonal = numpy.ones(n_hilbert)
linear_operator_diagonal = zeros_diagonal
# Loop through the terms.
for qubit_term in qubit_operator.terms:
is_zero = False
tensor_factor = 0
vecs = [ones_diagonal]
for pauli_operator in qubit_term:
op = pauli_operator[1]
if op in ['X', 'Y']:
is_zero = True
break
# Split vector by half and half for each bit.
if pauli_operator[0] > tensor_factor:
vecs = [v for iter_v in vecs for v in numpy.split(
iter_v, 2 ** (pauli_operator[0] - tensor_factor))]
vec_pairs = [numpy.split(v, 2) for v in vecs]
vecs = [v for vp in vec_pairs for v in (vp[0], -vp[1])]
tensor_factor = pauli_operator[0] + 1
if not is_zero:
linear_operator_diagonal += (qubit_operator.terms[qubit_term] *
numpy.concatenate(vecs))
return linear_operator_diagonal
def jw_configuration_state(occupied_orbitals, n_qubits):
"""Function to produce a basis state in the occupation number basis.
Args:
occupied_orbitals(list): A list of integers representing the indices
of the occupied orbitals in the desired basis state
n_qubits(int): The total number of qubits
Returns:
basis_vector(sparse): The basis state as a sparse matrix
"""
one_index = sum(2 ** (n_qubits - 1 - i) for i in occupied_orbitals)
basis_vector = numpy.zeros(2 ** n_qubits, dtype=float)
basis_vector[one_index] = 1
return basis_vector
def jw_hartree_fock_state(n_electrons, n_orbitals):
"""Function to produce Hartree-Fock state in JW representation."""
hartree_fock_state = jw_configuration_state(range(n_electrons),
n_orbitals)
return hartree_fock_state
def jw_number_indices(n_electrons, n_qubits):
"""Return the indices for n_electrons in n_qubits under JW encoding
Calculates the indices for all possible arrangements of n-electrons
within n-qubit orbitals when a Jordan-Wigner encoding is used.
Useful for restricting generic operators or vectors to a particular
particle number space when desired
Args:
n_electrons(int): Number of particles to restrict the operator to
n_qubits(int): Number of qubits defining the total state
Returns:
indices(list): List of indices in a 2^n length array that indicate
the indices of constant particle number within n_qubits
in a Jordan-Wigner encoding.
"""
occupations = itertools.combinations(range(n_qubits), n_electrons)
indices = [sum([2 ** n for n in occupation])
for occupation in occupations]
return indices
def jw_sz_indices(sz_value, n_qubits, n_electrons=None,
up_index=up_index, down_index=down_index):
r"""Return the indices of basis vectors with fixed Sz under JW encoding.
The returned indices label computational basis vectors which lie within
the corresponding eigenspace of the Sz operator,
.. math::
\begin{align}
S^{z} = \frac{1}{2}\sum_{i = 1}^{n}(n_{i, \alpha} - n_{i, \beta})
\end{align}
Args:
sz_value(float): Desired Sz value. Should be an integer or
half-integer.
n_qubits(int): Number of qubits defining the total state
n_electrons(int, optional): Number of particles to restrict the
operator to, if such a restriction is desired
up_index (Callable, optional): Function that maps a spatial index
to the index of the corresponding up site
down_index (Callable, optional): Function that maps a spatial index
to the index of the corresponding down site
Returns:
indices(list): The list of indices
"""
if n_qubits % 2 != 0:
raise ValueError('Number of qubits must be even')
if not (2. * sz_value).is_integer():
raise ValueError('Sz value must be an integer or half-integer')
n_sites = n_qubits // 2
sz_integer = int(2. * sz_value)
indices = []
if n_electrons is not None:
# Particle number is fixed, so the number of spin-up electrons
# (as well as the number of spin-down electrons) is fixed
if ((n_electrons + sz_integer) % 2 != 0 or
n_electrons < abs(sz_integer)):
raise ValueError('The specified particle number and sz value are '
'incompatible.')
num_up = (n_electrons + sz_integer) // 2
num_down = n_electrons - num_up
up_occupations = itertools.combinations(range(n_sites), num_up)
down_occupations = list(
itertools.combinations(range(n_sites), num_down))
# Each arrangement of up spins can be paired with an arrangement
# of down spins
for up_occupation in up_occupations:
up_occupation = [up_index(index) for index in up_occupation]
for down_occupation in down_occupations:
down_occupation = [down_index(index)
for index in down_occupation]
occupation = up_occupation + down_occupation
indices.append(sum(2 ** (n_qubits - 1 - k)
for k in occupation))
else:
# Particle number is not fixed
if sz_integer < 0:
# There are more down spins than up spins
more_map = down_index
less_map = up_index
else:
# There are at least as many up spins as down spins
more_map = up_index
less_map = down_index
for n in range(abs(sz_integer), n_sites + 1):
# Choose n of the 'more' spin and n - abs(sz_integer) of the
# 'less' spin
more_occupations = itertools.combinations(range(n_sites), n)
less_occupations = list(itertools.combinations(
range(n_sites), n - abs(sz_integer)))
# Each arrangement of the 'more' spins can be paired with an
# arrangement of the 'less' spin
for more_occupation in more_occupations:
more_occupation = [more_map(index)
for index in more_occupation]
for less_occupation in less_occupations:
less_occupation = [less_map(index)
for index in less_occupation]
occupation = more_occupation + less_occupation
indices.append(sum(2 ** (n_qubits - 1 - k)
for k in occupation))
return indices
def jw_number_restrict_operator(operator, n_electrons, n_qubits=None):
"""Restrict a Jordan-Wigner encoded operator to a given particle number
Args:
sparse_operator(ndarray or sparse): Numpy operator acting on
the space of n_qubits.
n_electrons(int): Number of particles to restrict the operator to
n_qubits(int): Number of qubits defining the total state
Returns:
new_operator(ndarray or sparse): Numpy operator restricted to
acting on states with the same particle number.
"""
if n_qubits is None:
n_qubits = int(numpy.log2(operator.shape[0]))
select_indices = jw_number_indices(n_electrons, n_qubits)
return operator[numpy.ix_(select_indices, select_indices)]
def jw_sz_restrict_operator(operator, sz_value,
n_electrons=None,
n_qubits=None,
up_index=up_index,
down_index=down_index):
"""Restrict a Jordan-Wigner encoded operator to a given Sz value
Args:
operator(ndarray or sparse): Numpy operator acting on
the space of n_qubits.
sz_value(float): Desired Sz value. Should be an integer or
half-integer.
n_electrons(int, optional): Number of particles to restrict the
operator to, if such a restriction is desired.
n_qubits(int, optional): Number of qubits defining the total state
up_index (Callable, optional): Function that maps a spatial index
to the index of the corresponding up site
down_index (Callable, optional): Function that maps a spatial index
to the index of the corresponding down site
Returns:
new_operator(ndarray or sparse): Numpy operator restricted to
acting on states with the desired Sz value.
"""
if n_qubits is None:
n_qubits = int(numpy.log2(operator.shape[0]))
select_indices = jw_sz_indices(
sz_value, n_qubits,
n_electrons=n_electrons,
up_index=up_index,
down_index=down_index)
return operator[numpy.ix_(select_indices, select_indices)]
def jw_number_restrict_state(state, n_electrons, n_qubits=None):
"""Restrict a Jordan-Wigner encoded state to a given particle number
Args:
state(ndarray or sparse): Numpy vector in
the space of n_qubits.
n_electrons(int): Number of particles to restrict the state to
n_qubits(int): Number of qubits defining the total state
Returns:
new_operator(ndarray or sparse): Numpy vector restricted to
states with the same particle number. May not be normalized.
"""
if n_qubits is None:
n_qubits = int(numpy.log2(state.shape[0]))
select_indices = jw_number_indices(n_electrons, n_qubits)
return state[select_indices]
def jw_sz_restrict_state(state, sz_value,
n_electrons=None,
n_qubits=None,
up_index=up_index,
down_index=down_index):
"""Restrict a Jordan-Wigner encoded state to a given Sz value
Args:
state(ndarray or sparse): Numpy vector in
the space of n_qubits.
sz_value(float): Desired Sz value. Should be an integer or
half-integer.
n_electrons(int, optional): Number of particles to restrict the
operator to, if such a restriction is desired.
n_qubits(int, optional): Number of qubits defining the total state
up_index (Callable, optional): Function that maps a spatial index
to the index of the corresponding up site
down_index (Callable, optional): Function that maps a spatial index
to the index of the corresponding down site
Returns:
new_operator(ndarray or sparse): Numpy vector restricted to
states with the desired Sz value. May not be normalized.
"""
if n_qubits is None:
n_qubits = int(numpy.log2(state.shape[0]))
select_indices = jw_sz_indices(
sz_value, n_qubits,
n_electrons=n_electrons,
up_index=up_index,
down_index=down_index)
return state[select_indices]
def jw_get_ground_state_at_particle_number(sparse_operator, particle_number):
"""Compute ground energy and state at a specified particle number.
Assumes the Jordan-Wigner transform. The input operator should be Hermitian
and particle-number-conserving.
Args:
sparse_operator(sparse): A Jordan-Wigner encoded sparse matrix.
particle_number(int): The particle number at which to compute the ground
energy and states
Returns:
ground_energy(float): The lowest eigenvalue of sparse_operator within
the eigenspace of the number operator corresponding to
particle_number.
ground_state(ndarray): The ground state at the particle number
"""
n_qubits = int(numpy.log2(sparse_operator.shape[0]))
# Get the operator restricted to the subspace of the desired particle number
restricted_operator = jw_number_restrict_operator(sparse_operator,
particle_number,
n_qubits)
# Compute eigenvalues and eigenvectors
if restricted_operator.shape[0] - 1 <= 1:
# Restricted operator too small for sparse eigensolver
dense_restricted_operator = restricted_operator.toarray()
eigvals, eigvecs = numpy.linalg.eigh(dense_restricted_operator)
else:
eigvals, eigvecs = scipy.sparse.linalg.eigsh(restricted_operator,
k=1,
which='SA')
# Expand the state
state = eigvecs[:, 0]
expanded_state = numpy.zeros(2 ** n_qubits, dtype=complex)
expanded_state[jw_number_indices(particle_number, n_qubits)] = state
return eigvals[0], expanded_state
def jw_get_gaussian_state(quadratic_hamiltonian, occupied_orbitals=None):
"""Compute an eigenvalue and eigenstate of a quadratic Hamiltonian.
Eigenstates of a quadratic Hamiltonian are also known as fermionic
Gaussian states.
Args:
quadratic_hamiltonian(QuadraticHamiltonian):
The Hamiltonian whose eigenstate is desired.
occupied_orbitals(list):
A list of integers representing the indices of the occupied
orbitals in the desired Gaussian state. If this is None
(the default), then it is assumed that the ground state is
desired, i.e., the orbitals with negative energies are filled.
Returns
-------
energy (float):
The eigenvalue.
state (sparse):
The eigenstate in scipy.sparse csc format.
"""
if not isinstance(quadratic_hamiltonian, QuadraticHamiltonian):
raise ValueError('Input must be an instance of QuadraticHamiltonian.')
n_qubits = quadratic_hamiltonian.n_qubits
# Compute the energy
orbital_energies, constant = quadratic_hamiltonian.orbital_energies()
if occupied_orbitals is None:
# The ground energy is desired
if quadratic_hamiltonian.conserves_particle_number:
num_negative_energies = numpy.count_nonzero(
orbital_energies < -EQ_TOLERANCE)
occupied_orbitals = range(num_negative_energies)
else:
occupied_orbitals = []
energy = numpy.sum(orbital_energies[occupied_orbitals]) + constant
# Obtain the circuit that prepares the Gaussian state
circuit_description, start_orbitals = gaussian_state_preparation_circuit(
quadratic_hamiltonian, occupied_orbitals)
# Initialize the starting state
state = jw_configuration_state(start_orbitals, n_qubits)
# Apply the circuit
if not quadratic_hamiltonian.conserves_particle_number:
particle_hole_transformation = (
jw_sparse_particle_hole_transformation_last_mode(n_qubits))
for parallel_ops in circuit_description:
for op in parallel_ops:
if op == 'pht':
state = particle_hole_transformation.dot(state)
else:
i, j, theta, phi = op
state = jw_sparse_givens_rotation(
i, j, theta, phi, n_qubits).dot(state)
return energy, state
def jw_slater_determinant(slater_determinant_matrix):
r"""Obtain a Slater determinant.
The input is an :math:`N_f \times N` matrix :math:`Q` with orthonormal
rows. Such a matrix describes the Slater determinant
.. math::
b^\dagger_1 \cdots b^\dagger_{N_f} \lvert \text{vac} \rangle,
where
.. math::
b^\dagger_j = \sum_{k = 1}^N Q_{jk} a^\dagger_k.
Args:
slater_determinant_matrix: The matrix :math:`Q` which describes the
Slater determinant to be prepared.
Returns:
The Slater determinant as a sparse matrix.
"""
circuit_description = slater_determinant_preparation_circuit(
slater_determinant_matrix)
start_orbitals = range(slater_determinant_matrix.shape[0])
n_qubits = slater_determinant_matrix.shape[1]
# Initialize the starting state
state = jw_configuration_state(start_orbitals, n_qubits)
# Apply the circuit
for parallel_ops in circuit_description:
for op in parallel_ops:
i, j, theta, phi = op
state = jw_sparse_givens_rotation(
i, j, theta, phi, n_qubits).dot(state)
return state
def jw_sparse_givens_rotation(i, j, theta, phi, n_qubits):
"""Return the matrix (acting on a full wavefunction) that performs a
Givens rotation of modes i and j in the Jordan-Wigner encoding."""
if j != i + 1:
raise ValueError('Only adjacent modes can be rotated.')
if j > n_qubits - 1:
raise ValueError('Too few qubits requested.')
cosine = numpy.cos(theta)
sine = numpy.sin(theta)
phase = numpy.exp(1.j * phi)
# Create the two-qubit rotation matrix
rotation_matrix = scipy.sparse.csc_matrix(
([1., phase * cosine, -phase * sine, sine, cosine, phase],
((0, 1, 1, 2, 2, 3), (0, 1, 2, 1, 2, 3))),
shape=(4, 4))
# Initialize identity operators
left_eye = scipy.sparse.eye(2 ** i, format='csc')
right_eye = scipy.sparse.eye(2 ** (n_qubits - 1 - j), format='csc')
# Construct the matrix and return
givens_matrix = kronecker_operators([left_eye, rotation_matrix, right_eye])
return givens_matrix
def jw_sparse_particle_hole_transformation_last_mode(n_qubits):
"""Return the matrix (acting on a full wavefunction) that performs a
particle-hole transformation on the last mode in the Jordan-Wigner
encoding.
"""
left_eye = scipy.sparse.eye(2 ** (n_qubits - 1), format='csc')
return kronecker_operators([left_eye, pauli_matrix_map['X']])
def get_density_matrix(states, probabilities):
n_qubits = states[0].shape[0]
density_matrix = scipy.sparse.csc_matrix(
(n_qubits, n_qubits), dtype=complex)
for state, probability in zip(states, probabilities):
state = scipy.sparse.csc_matrix(state.reshape((len(state), 1)))
density_matrix = density_matrix + probability * state * state.getH()
return density_matrix
def get_ground_state(sparse_operator, initial_guess=None):
"""Compute lowest eigenvalue and eigenstate.
Args:
sparse_operator (LinearOperator): Operator to find the ground state of.
initial_guess (ndarray): Initial guess for ground state. A good
guess dramatically reduces the cost required to converge.
Returns
-------
eigenvalue:
The lowest eigenvalue, a float.
eigenstate:
The lowest eigenstate in scipy.sparse csc format.
"""
values, vectors = scipy.sparse.linalg.eigsh(
sparse_operator, k=1, v0=initial_guess, which='SA', maxiter=1e7)
order = numpy.argsort(values)
values = values[order]
vectors = vectors[:, order]
eigenvalue = values[0]
eigenstate = vectors[:, 0]
return eigenvalue, eigenstate.T
def sparse_eigenspectrum(sparse_operator):
"""Perform a dense diagonalization.
Returns:
eigenspectrum: The lowest eigenvalues in a numpy array.
"""
dense_operator = sparse_operator.todense()
if is_hermitian(sparse_operator):
eigenspectrum = numpy.linalg.eigvalsh(dense_operator)
else:
eigenspectrum = numpy.linalg.eigvals(dense_operator)
return numpy.sort(eigenspectrum)
def expectation(operator, state):
"""Compute the expectation value of an operator with a state.
Args:
operator(scipy.sparse.spmatrix or scipy.sparse.linalg.LinearOperator):
The operator whose expectation value is desired.
state(numpy.ndarray or scipy.sparse.spmatrix): A numpy array
representing a pure state or a sparse matrix representing a density
matrix. If `operator` is a LinearOperator, then this must be a
numpy array.
Returns:
A complex number giving the expectation value.
Raises:
ValueError: Input state has invalid format.
"""
if isinstance(state, scipy.sparse.spmatrix):
# Handle density matrix.
if isinstance(operator, scipy.sparse.linalg.LinearOperator):
raise ValueError('Taking the expectation of a LinearOperator with '
'a density matrix is not supported.')
product = state * operator
expectation = numpy.sum(product.diagonal())
elif isinstance(state, numpy.ndarray):
# Handle state vector.
if len(state.shape) == 1:
# Row vector
expectation = numpy.dot(numpy.conjugate(state), operator * state)
else:
# Column vector
expectation = numpy.dot(numpy.conjugate(state.T),
operator * state)[0, 0]
else:
# Handle exception.
raise ValueError(
'Input state must be a numpy array or a sparse matrix.')
# Return.
return expectation
def variance(operator, state):
"""Compute variance of operator with a state.
Args:
operator(scipy.sparse.spmatrix or scipy.sparse.linalg.LinearOperator):
The operator whose expectation value is desired.
state(numpy.ndarray or scipy.sparse.spmatrix): A numpy array
representing a pure state or a sparse matrix representing a density
matrix.
Returns:
A complex number giving the variance.
Raises:
ValueError: Input state has invalid format.
"""
return (expectation(operator ** 2, state) -
expectation(operator, state) ** 2)
def expectation_computational_basis_state(operator, computational_basis_state):
"""Compute expectation value of operator with a state.
Args:
operator: Qubit or FermionOperator to evaluate expectation value of.
If operator is a FermionOperator, it must be normal-ordered.
computational_basis_state (scipy.sparse vector / list): normalized
computational basis state (if scipy.sparse vector), or list of
occupied orbitals.
Returns:
A real float giving expectation value.
Raises:
TypeError: Incorrect operator or state type.
"""
if isinstance(operator, QubitOperator):
raise NotImplementedError('Not yet implemented for QubitOperators.')
if not isinstance(operator, FermionOperator):
raise TypeError('operator must be a FermionOperator.')
occupied_orbitals = computational_basis_state
if not isinstance(occupied_orbitals, list):
computational_basis_state_index = (
occupied_orbitals.nonzero()[0][0])
occupied_orbitals = [digit == '1' for digit in
bin(computational_basis_state_index)[2:]][::-1]
expectation_value = operator.terms.get((), 0.0)
for i in range(len(occupied_orbitals)):
if occupied_orbitals[i]:
expectation_value += operator.terms.get(
((i, 1), (i, 0)), 0.0)
for j in range(i + 1, len(occupied_orbitals)):
expectation_value -= operator.terms.get(
((j, 1), (i, 1), (j, 0), (i, 0)), 0.0)
return expectation_value
def expectation_db_operator_with_pw_basis_state(
operator, plane_wave_occ_orbitals, n_spatial_orbitals, grid,
spinless):
"""Compute expectation value of a dual basis operator with a plane
wave computational basis state.
Args:
operator: Dual-basis representation of FermionOperator to evaluate
expectation value of. Can have at most 3-body terms.
plane_wave_occ_orbitals (list): list of occupied plane-wave orbitals.
n_spatial_orbitals (int): Number of spatial orbitals.
grid (openfermion.utils.Grid): The grid used for discretization.
spinless (bool): Whether the system is spinless.
Returns:
A real float giving the expectation value.
"""
expectation_value = operator.terms.get((), 0.0)
for single_action, coefficient in operator.terms.items():
if len(single_action) == 2:
expectation_value += coefficient * (
expectation_one_body_db_operator_computational_basis_state(
single_action, plane_wave_occ_orbitals, grid, spinless) /
n_spatial_orbitals)
elif len(single_action) == 4:
expectation_value += coefficient * (
expectation_two_body_db_operator_computational_basis_state(
single_action, plane_wave_occ_orbitals, grid, spinless) /
n_spatial_orbitals ** 2)
elif len(single_action) == 6:
expectation_value += coefficient * (
expectation_three_body_db_operator_computational_basis_state(
single_action, plane_wave_occ_orbitals, grid, spinless) /
n_spatial_orbitals ** 3)
return expectation_value
def expectation_one_body_db_operator_computational_basis_state(
dual_basis_action, plane_wave_occ_orbitals, grid, spinless):
"""Compute expectation value of a 1-body dual-basis operator with a
plane wave computational basis state.
Args:
dual_basis_action: Dual-basis action of FermionOperator to
evaluate expectation value of.
plane_wave_occ_orbitals (list): list of occupied plane-wave orbitals.
grid (openfermion.utils.Grid): The grid used for discretization.
spinless (bool): Whether the system is spinless.
Returns:
A real float giving the expectation value.
"""
expectation_value = 0.0
r_p = grid.position_vector(grid.grid_indices(dual_basis_action[0][0],
spinless))
r_q = grid.position_vector(grid.grid_indices(dual_basis_action[1][0],
spinless))
for orbital in plane_wave_occ_orbitals:
# If there's spin, p and q have to have the same parity (spin),
# and the new orbital has to have the same spin as these.
k_orbital = grid.momentum_vector(grid.grid_indices(orbital,
spinless))
# The Fourier transform is spin-conserving. This means that p, q,
# and the new orbital all have to have the same spin (parity).
if spinless or (dual_basis_action[0][0] % 2 ==
dual_basis_action[1][0] % 2 == orbital % 2):
expectation_value += numpy.exp(-1j * k_orbital.dot(r_p - r_q))
return expectation_value
def expectation_two_body_db_operator_computational_basis_state(
dual_basis_action, plane_wave_occ_orbitals, grid, spinless):
"""Compute expectation value of a 2-body dual-basis operator with a
plane wave computational basis state.
Args:
dual_basis_action: Dual-basis action of FermionOperator to
evaluate expectation value of.
plane_wave_occ_orbitals (list): list of occupied plane-wave orbitals.
grid (openfermion.utils.Grid): The grid used for discretization.
spinless (bool): Whether the system is spinless.
Returns:
A float giving the expectation value.
"""
expectation_value = 0.0
r = {}
for i in range(4):
r[i] = grid.position_vector(grid.grid_indices(dual_basis_action[i][0],
spinless))
rr = {}
k_map = {}
for i in range(2):
rr[i] = {}
k_map[i] = {}
for j in range(2, 4):
rr[i][j] = r[i] - r[j]
k_map[i][j] = {}
# Pre-computations.
for o in plane_wave_occ_orbitals:
k = grid.momentum_vector(grid.grid_indices(o, spinless))
for i in range(2):
for j in range(2, 4):
k_map[i][j][o] = k.dot(rr[i][j])
for orbital1 in plane_wave_occ_orbitals:
k1ac = k_map[0][2][orbital1]
k1ad = k_map[0][3][orbital1]
for orbital2 in plane_wave_occ_orbitals:
if orbital1 != orbital2:
k2bc = k_map[1][2][orbital2]
k2bd = k_map[1][3][orbital2]
# The Fourier transform is spin-conserving. This means that
# the parity of the orbitals involved in the transition must
# be the same.
if spinless or (
(dual_basis_action[0][0] % 2 ==
dual_basis_action[3][0] % 2 == orbital1 % 2) and
(dual_basis_action[1][0] % 2 ==
dual_basis_action[2][0] % 2 == orbital2 % 2)):
value = numpy.exp(-1j * (k1ad + k2bc))
# Add because it came from two anti-commutations.
expectation_value += value
# The Fourier transform is spin-conserving. This means that
# the parity of the orbitals involved in the transition must
# be the same.
if spinless or (
(dual_basis_action[0][0] % 2 ==
dual_basis_action[2][0] % 2 == orbital1 % 2) and
(dual_basis_action[1][0] % 2 ==
dual_basis_action[3][0] % 2 == orbital2 % 2)):
value = numpy.exp(-1j * (k1ac + k2bd))
# Subtract because it came from a single anti-commutation.
expectation_value -= value
return expectation_value
def expectation_three_body_db_operator_computational_basis_state(
dual_basis_action, plane_wave_occ_orbitals, grid, spinless):
"""Compute expectation value of a 3-body dual-basis operator with a
plane wave computational basis state.
Args:
dual_basis_action: Dual-basis action of FermionOperator to
evaluate expectation value of.
plane_wave_occ_orbitals (list): list of occupied plane-wave orbitals.
grid (openfermion.utils.Grid): The grid used for discretization.
spinless (bool): Whether the system is spinless.
Returns:
A float giving the expectation value.
"""
expectation_value = 0.0
r = {}
for i in range(6):
r[i] = grid.position_vector(grid.grid_indices(dual_basis_action[i][0],
spinless))
rr = {}
k_map = {}
for i in range(3):
rr[i] = {}
k_map[i] = {}
for j in range(3, 6):
rr[i][j] = r[i] - r[j]
k_map[i][j] = {}
# Pre-computations.
for o in plane_wave_occ_orbitals:
k = grid.momentum_vector(grid.grid_indices(o, spinless))
for i in range(3):
for j in range(3, 6):
k_map[i][j][o] = k.dot(rr[i][j])
for orbital1 in plane_wave_occ_orbitals:
k1ad = k_map[0][3][orbital1]
k1ae = k_map[0][4][orbital1]
k1af = k_map[0][5][orbital1]
for orbital2 in plane_wave_occ_orbitals:
if orbital1 != orbital2:
k2bd = k_map[1][3][orbital2]
k2be = k_map[1][4][orbital2]
k2bf = k_map[1][5][orbital2]
for orbital3 in plane_wave_occ_orbitals:
if orbital1 != orbital3 and orbital2 != orbital3:
k3cd = k_map[2][3][orbital3]
k3ce = k_map[2][4][orbital3]
k3cf = k_map[2][5][orbital3]
# Handle \delta_{ad} \delta_{bf} \delta_{ce} after FT.
# The Fourier transform is spin-conserving.
if spinless or (
(dual_basis_action[0][0] % 2 ==
dual_basis_action[3][0] % 2 ==
orbital1 % 2) and
(dual_basis_action[1][0] % 2 ==
dual_basis_action[5][0] % 2 ==
orbital2 % 2) and
(dual_basis_action[2][0] % 2 ==
dual_basis_action[4][0] % 2 ==
orbital3 % 2)):
expectation_value += numpy.exp(-1j * (
k1ad + k2bf + k3ce))
# Handle -\delta_{ad} \delta_{be} \delta_{cf} after FT.
# The Fourier transform is spin-conserving.
if spinless or (
(dual_basis_action[0][0] % 2 ==
dual_basis_action[3][0] % 2 ==
orbital1 % 2) and
(dual_basis_action[1][0] % 2 ==
dual_basis_action[4][0] % 2 ==
orbital2 % 2) and
(dual_basis_action[2][0] % 2 ==
dual_basis_action[5][0] % 2 ==
orbital3 % 2)):
expectation_value -= numpy.exp(-1j * (
k1ad + k2be + k3cf))
# Handle -\delta_{ae} \delta_{bf} \delta_{cd} after FT.
# The Fourier transform is spin-conserving.
if spinless or (
(dual_basis_action[0][0] % 2 ==
dual_basis_action[4][0] % 2 ==
orbital1 % 2) and
(dual_basis_action[1][0] % 2 ==
dual_basis_action[5][0] % 2 ==
orbital2 % 2) and
(dual_basis_action[2][0] % 2 ==
dual_basis_action[3][0] % 2 ==
orbital3 % 2)):
expectation_value -= numpy.exp(-1j * (
k1ae + k2bf + k3cd))
# Handle \delta_{ae} \delta_{bd} \delta_{cf} after FT.
# The Fourier transform is spin-conserving.
if spinless or (
(dual_basis_action[0][0] % 2 ==
dual_basis_action[4][0] % 2 ==
orbital1 % 2) and
(dual_basis_action[1][0] % 2 ==
dual_basis_action[3][0] % 2 ==
orbital2 % 2) and
(dual_basis_action[2][0] % 2 ==
dual_basis_action[5][0] % 2 ==
orbital3 % 2)):
expectation_value += numpy.exp(-1j * (
k1ae + k2bd + k3cf))
# Handle \delta_{af} \delta_{be} \delta_{cd} after FT.
# The Fourier transform is spin-conserving.
if spinless or (
(dual_basis_action[0][0] % 2 ==
dual_basis_action[5][0] % 2 ==
orbital1 % 2) and
(dual_basis_action[1][0] % 2 ==
dual_basis_action[4][0] % 2 ==
orbital2 % 2) and
(dual_basis_action[2][0] % 2 ==
dual_basis_action[3][0] % 2 ==
orbital3 % 2)):
expectation_value += numpy.exp(-1j * (
k1af + k2be + k3cd))
# Handle -\delta_{af} \delta_{bd} \delta_{ce} after FT.
# The Fourier transform is spin-conserving.
if spinless or (
(dual_basis_action[0][0] % 2 ==
dual_basis_action[5][0] % 2 ==
orbital1 % 2) and
(dual_basis_action[1][0] % 2 ==
dual_basis_action[3][0] % 2 ==
orbital2 % 2) and
(dual_basis_action[2][0] % 2 ==
dual_basis_action[4][0] % 2 ==
orbital3 % 2)):
expectation_value -= numpy.exp(-1j * (
k1af + k2bd + k3ce))
return expectation_value
def get_gap(sparse_operator, initial_guess=None):
"""Compute gap between lowest eigenvalue and first excited state.
Args:
sparse_operator (LinearOperator): Operator to find the ground state of.
initial_guess (ndarray): Initial guess for eigenspace. A good
guess dramatically reduces the cost required to converge.
Returns: A real float giving eigenvalue gap.
"""
if not is_hermitian(sparse_operator):
raise ValueError('sparse_operator must be Hermitian.')
values, _ = scipy.sparse.linalg.eigsh(
sparse_operator, k=2, v0=initial_guess, which='SA', maxiter=1e7)
gap = abs(values[1] - values[0])
return gap
def inner_product(state_1, state_2):
"""Compute inner product of two states."""
return numpy.dot(state_1.conjugate(), state_2)
def boson_ladder_sparse(n_modes, mode, ladder_type, trunc):
r"""Make a matrix representation of a singular bosonic ladder operator
in the Fock space.
Since the bosonic operator lies in an infinite Fock space,
a truncation value needs to be provide so that a sparse matrix
of finite size can be returned.
Args:
n_modes (int): Number of modes in the system Hilbert space.
mode (int): The mode the ladder operator targets.
ladder_type (int): This is a nonzero integer. 0 indicates a lowering
operator, 1 a raising operator.
trunc (int): The size at which the Fock space should be truncated
when returning the matrix representing the ladder operator.
Returns:
The corresponding trunc x trunc Scipy sparse matrix.
"""
if trunc < 1 or not isinstance(trunc, int):
raise ValueError("Fock space truncation must be a positive integer.")
if ladder_type:
lop = scipy.sparse.spdiags(numpy.sqrt(range(1, trunc)),
-1, trunc, trunc, format='csc')
else:
lop = scipy.sparse.spdiags(numpy.sqrt(range(trunc)),
1, trunc, trunc, format='csc')
Id = [scipy.sparse.identity(trunc, format='csc', dtype=complex)]
operator_list = Id*mode + [lop] + Id*(n_modes - mode - 1)
operator = kronecker_operators(operator_list)
return operator
def single_quad_op_sparse(n_modes, mode, quadrature, hbar, trunc):
r"""Make a matrix representation of a singular quadrature
operator in the Fock space.
Since the bosonic operators lie in an infinite Fock space,
a truncation value needs to be provide so that a sparse matrix
of finite size can be returned.
Args:
n_modes (int): Number of modes in the system Hilbert space.
mode (int): The mode the ladder operator targets.
quadrature (str): 'q' for the canonical position operator,
'p' for the canonical moment]um operator.
hbar (float): the value of hbar to use in the definition of the
canonical commutation relation [q_i, p_j] = \delta_{ij} i hbar.
trunc (int): The size at which the Fock space should be truncated
when returning the matrix representing the ladder operator.
Returns:
The corresponding trunc x trunc Scipy sparse matrix.
"""
if trunc < 1 or not isinstance(trunc, int):
raise ValueError("Fock space truncation must be a positive integer.")
b = boson_ladder_sparse(1, 0, 0, trunc)
if quadrature == 'q':
op = numpy.sqrt(hbar/2) * (b + b.conj().T)
elif quadrature == 'p':
op = -1j*numpy.sqrt(hbar/2) * (b - b.conj().T)
Id = [scipy.sparse.identity(trunc, dtype=complex, format='csc')]
operator_list = Id*mode + [op] + Id*(n_modes - mode - 1)
operator = kronecker_operators(operator_list)
return operator
def boson_operator_sparse(operator, trunc, hbar=1.):
r"""Initialize a Scipy sparse matrix in the Fock space
from a bosonic operator.
Since the bosonic operators lie in an infinite Fock space,
a truncation value needs to be provide so that a sparse matrix
of finite size can be returned.
Args:
operator: One of either BosonOperator or QuadOperator.
trunc (int): The size at which the Fock space should be truncated
when returning the matrix representing the ladder operator.
hbar (float): the value of hbar to use in the definition of the
canonical commutation relation [q_i, p_j] = \delta_{ij} i hbar.
This only applies if calcualating the sparse representation of
a quadrature operator.
Returns:
The corresponding Scipy sparse matrix of size [trunc, trunc].
"""
if isinstance(operator, QuadOperator):
from openfermion.transforms._conversion import get_boson_operator
boson_operator = get_boson_operator(operator, hbar)
elif isinstance(operator, BosonOperator):
boson_operator = operator
else:
raise ValueError("Only BosonOperator and QuadOperator are supported.")
if trunc < 1 or not isinstance(trunc, int):
raise ValueError("Fock space truncation must be a positive integer.")
# count the number of modes
n_modes = 0
for term in boson_operator.terms:
for ladder_operator in term:
if ladder_operator[0] + 1 > n_modes:
n_modes = ladder_operator[0] + 1
# Construct the Scipy sparse matrix.
n_hilbert = trunc ** n_modes
values_list = [[]]
row_list = [[]]
column_list = [[]]
# Loop through the terms.
for term in boson_operator.terms:
coefficient = boson_operator.terms[term]
term_operator = coefficient*scipy.sparse.identity(
n_hilbert, dtype=complex, format='csc')
for ladder_op in term:
# Add actual operator to the list.
b = boson_ladder_sparse(n_modes, ladder_op[0], ladder_op[1], trunc)
term_operator = term_operator.dot(b)
# Extract triplets from sparse_term.
values_list.append(term_operator.tocoo(copy=False).data)
(row, column) = term_operator.nonzero()
column_list.append(column)
row_list.append(row)
# Create sparse operator.
values_list = numpy.concatenate(values_list)
row_list = numpy.concatenate(row_list)
column_list = numpy.concatenate(column_list)
sparse_operator = scipy.sparse.coo_matrix((
values_list, (row_list, column_list)),
shape=(n_hilbert, n_hilbert)).tocsc(copy=False)
sparse_operator.eliminate_zeros()
return sparse_operator
