Subspace Hamiltonians

Once one has a qubit operator and subspace in hand, it is time to construct the corresponding SubspaceHamiltonian object that is used to build matrix representations of the Hamiltonian confied to the given subspace or provide the matrix-vector product operation for matrix-free evaluation.

import numpy as np
import scipy.sparse.linalg as spla
import fulqrum as fq

To begin, let us import a Fermionic operator for LiH from a JSON file and convert it to a qubit operator:

fop = fq.FermionicOperator.from_json("./data/lih.json")
op = fop.extended_jw_transformation()

This operator is small, consisting of only 12 qubits:

op.width

12

so we can construct the subspace that spans the full Hilbert space:

counts = []
for kk in range(2**op.width):
    counts.append(bin(kk)[2:].zfill(op.width))

S = fq.Subspace([counts])

Constucting subspace Hamiltonians

To build a SubspaceHamiltonian object one needs, at minimum, a qubit operator (Hamiltonian) and usually a corresponding subspace in which the operator is to be projected:

Hsub = fq.SubspaceHamiltonian(op, S)

when this object is constructed a collection of transformations are applied to the input Hamiltonian that optimize it for evaluating its matrix-elements efficiently. This object is also a subclass of the scipy.sparse.linalg.LinearOperator class, allowing it to be used for matrix-free computations in those routines that support this interface.

Using subspace Hamiltonians

Having built a SubspaceHamiltonian we can now put it to work to solve eigenvalue equations of interest, e.g. finding the gound state energy of our molecule.

To begin, lets construct a CSR matrix for the Hamiltonian within the subspace:

A = Hsub.to_csr_linearoperator()

Note that we call this a “linearoperator” because it is actually a scipy.sparse.linalg.LinearOperator wrapper around the CSR matrix.

A.matrix

<Compressed Sparse Row sparse array of dtype 'float64'
        with 102400 stored elements and shape (4096, 4096)>

This allows us to perform the matrix-vector product using OpenMP, and thus accelerate the operation. This can then be passed onto a classical eigensolver for evaluation

evals, evecs = spla.eigsh(A, k=1, which="SA")
evals

array([-7.87565256])

There is also a “fast” version of this CSR matrix building, that is faster than the standard method, at the expense of using 2x the memory. This is likely the best method to try first

B = Hsub.to_csr_linearoperator_fast()

evals, evecs = spla.eigsh(B, k=1, which="SA")
evals

array([-7.87565256])

Finally, we can pass the entire SubspaceHamiltonian to the eigensolver itself to perform matrix-free evaluation. This avoids any memory allocation for the operator itself, but comes at the cost of extended runtime. The result of course will be identical:

evals, evecs = spla.eigsh(Hsub, k=1, which="SA")
evals

array([-7.87565256])

Extracting eigenvectors

After one or more eigenvectors is found they are expressed in an ordering that is determined by the order of the bit-strings in the Subspace. As this ordering may not be apparent, it is best to use interpret_vector to get the dictionary representation of the eigenvectors

Hsub.interpret_vector(evecs, atol=1e-3)

{'000011000011': 0.991017311745354,
 '000011000101': 0.029187947808189257,
 '000011100100': 0.0026739729129043027,
 '000101000011': 0.029187947808187394,
 '000101000101': -0.01982607762509705,
 '000101100001': 0.042018939160017275,
 '000101100010': -0.0012816405415078255,
 '000110000110': -0.0035432760522673894,
 '000110100001': -0.003955613454412125,
 '000110100010': -0.0014842102082854806,
 '001001001001': -0.035699337389609336,
 '001001001010': 0.002926474405920963,
 '001010001001': 0.002926474405920973,
 '001010001010': -0.001741011876849469,
 '010001010001': -0.035699337389608496,
 '010001010010': 0.002926474405920986,
 '010010010001': 0.0029264744059209827,
 '010010010010': -0.001741011876849485,
 '100001000101': 0.04201893916001825,
 '100001000110': -0.00395561345441214,
 '100001100001': -0.0979081466961776,
 '100001100010': 0.002784758898105246,
 '100010000101': -0.0012816405415078276,
 '100010000110': -0.0014842102082855213,
 '100010100001': 0.00278475889810527,
 '100100000011': 0.0026739729129043053}

Here we use the atol keyword argument to select only those bit-strings with amplitudes greater than the absolute tolerance value. Recall that the absolute value of the amplitudes squared is the corresponding probability.

Updating the subspace in a SubspaceHamiltonian

Initialization of a SubspaceHamiltonian performs several transformations on the input Hamiltonian. These transformations take some amount of time, and this time can add up if using multiple subspaces for the same operator. To avoid this overhead, it is best to update the subspace for a pre-existing SubspaceHamiltonian, as opposed to creating an entire new one.

Tol show how this works, lets create a subspace Hamiltonian with no subspace:

Hsub2 = fq.SubspaceHamiltonian(op)

It is impossible to use Hsub2 as it stands now, but we can update the subspace and use it like before:

Hsub2.update_subspace(S)

evals, evecs = spla.eigsh(Hsub2, k=1, which="SA")
evals

array([-7.87565256])

The subspace can be updated as many times as needed.