Reduced Density Matrices

Reduced Density Matrices (RDMs) are a use case where measurement reduction in combination with computables is very useful. Naively, we can evaluate each term in the n-body RDM by evaluating many ExpectationValue computables independently on each spin-traced excitation operator. However, across the entire RDM, this would result in many redundant measurements, making it advantageous to create a single, composite computable for an RDM ndarray and use a protocol to build the optimal number of circuits for the whole RDM.

To build a 1-RDM operator, we can use InQuanto’s FermionSpace:

import numpy
from inquanto.spaces import FermionSpace

space = FermionSpace(4)
rdm_operators = space.construct_one_body_spatial_rdm_operators((FermionSpace.SPIN_UP, FermionSpace.SPIN_UP))
print(rdm_operators)
print(type(rdm_operators))

[[{((0, 1), (0, 0)): 1.0} {((0, 1), (2, 0)): 1.0}]
 [{((2, 1), (0, 0)): 1.0} {((2, 1), (2, 0)): 1.0}]]
<class 'numpy.ndarray'>

As a result, we obtain a numpy ndarray of FermionOperator objects. To compute the expectation values, we also need a specific state. In this example, we will use the TrotterAnsatz:

from inquanto.ansatzes import TrotterAnsatz
from inquanto.operators import QubitOperator, QubitOperatorList
from inquanto.states import QubitState

ansatz = TrotterAnsatz(
    QubitOperatorList.from_list([QubitOperator("Y0 X1 X2 X3", 1j)]),
    QubitState([1, 1, 0, 0]),
)

parameters = ansatz.state_symbols.construct_from_array([0.11])

With the ansatz and RDM operators, we can construct a single computable for the RDM:

from inquanto.computables import ExpectationValueNonHermitian
from inquanto.computables.primitive import ComputableNDArray

# We make a broadcast function over the ndarray that converts every element into
# an ExpectationValueNonHermitian computable
ndarray_broadcast = numpy.vectorize(
    lambda operator: ExpectationValueNonHermitian(ansatz, operator.qubit_encode())
)

# Then, we can convert the ndarray computables to a ComputableNDArray of computables
# to ensure that we can run and evaluate recursively.
rdm_computable = ComputableNDArray(ndarray_broadcast(rdm_operators))

print(rdm_computable)
print(type(rdm_computable))

[[ExpectationValueNonHermitian(state=<inquanto.ansatzes._trotter_ansatz.TrotterAnsatz, qubits=4, gates=3, symbols=1>, kernel={(): np.float64(0.5), (Zq[0]): np.float64(-0.5)})
  ExpectationValueNonHermitian(state=<inquanto.ansatzes._trotter_ansatz.TrotterAnsatz, qubits=4, gates=3, symbols=1>, kernel={(Yq[0], Zq[1], Xq[2]): (-0-0.25j), (Yq[0], Zq[1], Yq[2]): 0.25, (Xq[0], Zq[1], Xq[2]): 0.25, (Xq[0], Zq[1], Yq[2]): 0.25j})]
 [ExpectationValueNonHermitian(state=<inquanto.ansatzes._trotter_ansatz.TrotterAnsatz, qubits=4, gates=3, symbols=1>, kernel={(Yq[0], Zq[1], Xq[2]): 0.25j, (Xq[0], Zq[1], Xq[2]): 0.25, (Yq[0], Zq[1], Yq[2]): 0.25, (Xq[0], Zq[1], Yq[2]): (-0-0.25j)})
  ExpectationValueNonHermitian(state=<inquanto.ansatzes._trotter_ansatz.TrotterAnsatz, qubits=4, gates=3, symbols=1>, kernel={(): 0.5, (Zq[2]): -0.5})]]
<class 'inquanto.computables.primitive.collections._ndarray.ComputableNDArray'>

With this RDM computable, we can use protocols to create an evaluator function. We first use default_evaluate() to perform a statevector calculation as a test:

print(rdm_computable.default_evaluate(parameters))

[[(0.9879487246653031+0j) (1.7242389772158059e-18+3.489965711775505e-32j)]
 [(1.7242389772158059e-18-3.489965711775505e-32j)
  (0.012051275334697192+0j)]]

Now, we use the shot-based PauliAveraging protocol to evaluate the RDM computable:

from pytket.extensions.qiskit import AerBackend
from inquanto.protocols import PauliAveraging
from pytket.partition import PauliPartitionStrat

# We can use the PauliAveraging protocol to get an evaluator
evaluator = (
    PauliAveraging(
        AerBackend(),
        shots_per_circuit = 10000,
        pauli_partition_strategy=PauliPartitionStrat.CommutingSets
    )
    .build_from(parameters, rdm_computable)
    .compile_circuits()
    .run(seed=0)
    .get_evaluator()
)

rdm_measured = rdm_computable.evaluate(evaluator)

print(rdm_measured)

[[(0.9867+0j) (0.0029+0.0041j)]
 [(0.0029-0.0041j) (0.013299999999999979+0j)]]

InQuanto provides built-in support for RDM computables such as the above, for example with the RestrictedOneBodyRDMComputable class, which evaluates to a RestrictedOneBodyRDM. A suite of other density matrix composite computables are supported, for instance, the SpinlessNBodyRDMArrayRealComputable can build n-body RDMs and reduce measurement circuits further by taking symmetries into account:

from inquanto.ansatzes import FermionSpaceAnsatzChemicallyAwareUCCSD
from inquanto.computables.composite import SpinlessNBodyRDMArrayRealComputable
from inquanto.express import load_h5
from inquanto.mappings import QubitMappingJordanWigner

from inquanto.spaces import QubitSpace, FermionSpace
from inquanto.states import FermionState

qubit_hamiltonian = load_h5(
    "h2_631g_symmetry.h5", as_tuple=True
).hamiltonian_operator.qubit_encode()

space = FermionSpace(
    8,
    point_group="D2h",
    orb_irreps=numpy.asarray(["Ag", "Ag", "B1u", "B1u", "Ag", "Ag", "B1u", "B1u"]),
)

ansatz = FermionSpaceAnsatzChemicallyAwareUCCSD(
    fermion_space=space, fermion_state=FermionState([1, 1, 0, 0, 0, 0, 0, 0])
)

qubit_space = QubitSpace(space.n_spin_orb)

symmetry_operators = qubit_space.symmetry_operators_z2(qubit_hamiltonian)

rdm1 = SpinlessNBodyRDMArrayRealComputable(
    1,
    space,
    ansatz,
    QubitMappingJordanWigner(),
    symmetry_operators,
)

print(rdm1)

<inquanto.computables.composite.rdm._rdm_real.SpinlessNBodyRDMArrayRealComputable object at 0x7f4664d92850>

Above, we used the express module to load a Hamiltonian, and use a Fermion space to generate the Z2 symmetries of our Hamiltonian. The SpinlessNBodyRDMArrayRealComputable uses the symmetry information to filter for Z2 symmetry.