Overlap Matrices & Non-Orthogonal Subspaces

InQuanto provides the composite computable OverlapMatrixComputable for calculating overlap matrices:

(33)\[S_{ij} = \langle \Psi_i | \Psi_j \rangle\]

given a list of ansatzes \(\{|\Psi_0\rangle, |\Psi_1\rangle, \dots\}\). It also supports matrices of overlaps with a Hermitian kernel:

(34)\[O_{ij} = \langle \Psi_i | \hat{O} | \Psi_j \rangle\]

where \(\hat{O}\) is represented by a QubitOperator. Focusing on the latter of these objects, we first show a statevector calculation with a simple kernel:

from inquanto.ansatzes import TrotterAnsatz
from inquanto.states import QubitState
from inquanto.operators import QubitOperatorList, QubitOperator
from inquanto.computables.composite import OverlapMatrixComputable
from inquanto.core import SymbolDict

states = [
    TrotterAnsatz(
        QubitOperatorList.from_string("a [(1j, Y0 X1 X2 X3)]"),
        QubitState([1, 1, 0, 0])
    ),
    TrotterAnsatz(
        QubitOperatorList.from_string("b [(1j, Y0 Z1 Z2 Z3)]"),
        QubitState([1, 1, 0, 0])
    ),
]
kernel = QubitOperator.from_string("(1, Z0 Z1)")
params = SymbolDict(a=0.5, b=0.5)

om_computable = OverlapMatrixComputable(states, kernel)
om_result_statevector = om_computable.default_evaluate(params)
print(om_result_statevector)
[[1.  +0.j 0.77+0.j]
 [0.77-0.j 0.54+0.j]]

This computable uses ExpectationValue computables for the diagonal elements of the overlap matrix, and Overlap computables for the off-diagonal elements. Since the kernel must be Hermitian, it is true that \(O_{ij} = O_{ji}^\dagger\). The OverlapMatrixComputable uses this symmetry to this result to reduce the total number of calculations required.

To perform a shot-based experiment, we require two different InQuanto protocols, one for calculating the ExpectationValue components, and one for the Overlap components. In this example we choose the PauliAveraging and HadamardTestOverlap protocols for these respective tasks. See the Protocols manual page for more information.

Below we use the build_protocols_from() method to parse the overlap matrix computable and construct a ProtocolList object. The ProtocolList class groups shot-based protocols together to measure a composite computable:

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

shot_backend = AerBackend()

expval_protocols = PauliAveraging.build_protocols_from(
    parameters=params,
    computable=om_computable,
    backend=shot_backend,
    shots_per_circuit=10000,
    pauli_partition_strategy=PauliPartitionStrat.CommutingSets,
)

overlap_protocols = HadamardTestOverlap.build_protocols_from(
    parameters=params,
    computable=om_computable,
    backend=shot_backend,
    shots_per_circuit=10000,
    direct=True,
    pauli_partition_strategy=PauliPartitionStrat.CommutingSets,
)

protocol_list = expval_protocols + overlap_protocols
print(protocol_list.dataframe_protocol_circuit())
       Protocol ID        Protocol Type  Qubits  Depth2q  Shots
0  140450607585680       PauliAveraging       4        6  10000
1  140451815130960       PauliAveraging       4        6  10000
2  140450606514512  HadamardTestOverlap       5       87  10000
3  140450606514512  HadamardTestOverlap       5       87  10000

The dataframe_protocol_circuit() method shows a breakdown of the circuits required to measure the overlap matrix. Two circuits come from the PauliAveraging protocol for measuring the expectation values on the diagonal, and two more from the HadamardTestOverlap protocol for measuring the real and imaginary parts of the single, unique off-diagonal element.

Below, we run these circuits using the pytket AerBackend, evaluate the overlap matrix, and compare to the statevector result obtained above:

protocol_list.run(seed=0)
om_result_shotbased = om_computable.evaluate(protocol_list.get_evaluator())

print(f"Statevector result: \n{om_result_statevector}\n")
print(f"Shot-based result: \n{om_result_shotbased}")
Statevector result: 
[[1.  +0.j 0.77+0.j]
 [0.77-0.j 0.54+0.j]]

Shot-based result: 
[[1.   +0.j   0.772-0.01j]
 [0.772+0.01j 0.53 +0.j  ]]

InQuanto also provides the specialized composite computable: NonOrthogonalMatricesComputable for building the generalized eigenvalue problem in a non-orthogonal subspace:

(35)\[HC = SCE\]

where \(H\) is the hamiltonian matrix and \(S\) is the overlap matrix in some subspace \(\{|\Psi_i\rangle\}\):

(36)\[H_{ij} = \langle \Psi_i | \hat{H} | \Psi_j \rangle\]
(37)\[S_{ij} = \langle \Psi_i | \Psi_j \rangle.\]

Given a set of states, the NonOrthogonalMatricesComputable constructs both such matrices using OverlapMatrixComputables discussed above. Once evaluated, we may solve the eigenvalue problem straightforwardly with the pd_safe_eigh() function from inquanto.core. The ground state energy obtained in this approach is bounded from below by the configuration interaction (CI) energy.

Below we demonstrate a simple calculation of H2 in a minimal basis, using the two Trotter states defined above as a crude, non-orthogonal subspace:

from inquanto.computables.composite import NonOrthogonalMatricesComputable
from inquanto.express import load_h5
from inquanto.core import pd_safe_eigh

h2_data = load_h5("h2_sto3g.h5")
qubit_hamiltonian = h2_data["hamiltonian_operator"].qubit_encode()

no_computable = NonOrthogonalMatricesComputable(
    qubit_hamiltonian,
    states
)

# Using the same ansatz parameters as above
expval_protocols = PauliAveraging.build_protocols_from(
    parameters=params,
    computable=no_computable,
    backend=shot_backend,
    shots_per_circuit=10000,
    pauli_partition_strategy=PauliPartitionStrat.CommutingSets,
)

overlap_protocols = HadamardTestOverlap.build_protocols_from(
    parameters=params,
    computable=no_computable,
    backend=shot_backend,
    shots_per_circuit=10000,
    direct=True,
    pauli_partition_strategy=PauliPartitionStrat.CommutingSets,
)

no_protocols = expval_protocols + overlap_protocols
no_protocols.run(seed=0)

h, s = no_computable.evaluate(no_protocols.get_evaluator())
e, _, _ = pd_safe_eigh(h, s)

print(f"CASCI energy: {h2_data['energy_casci']}")
print(f"NO shot-based energy: {e[0]}")
CASCI energy: -1.1368465754720547
NO shot-based energy: -0.9825372067689355