SubspaceHamiltonian sometimes drops terms from extended JW operators

[willkirby1]

Summary

SubspaceHamiltonian drops a valid weight-4 off-diagonal term produced by FermionicOperator.extended_jw_transformation().

The direct QubitOperator.matrix_element() calculation is correct, but the matrix-free, CSR, and fast CSR SubspaceHamiltonian paths all produce a zero matrix.

Reproduced against current main at commit dafb4091d552d130c42d13dc74a3484e12d70d3f with Fulqrum 0.3.0.

Minimal reproducer

import numpy as np
import fulqrum as fq

fermionic_op = fq.FermionicOperator(
    28,
    [
        ("+-+-", [16, 0, 15, 1], 1.0),
        ("+-+-", [1, 15, 0, 16], 1.0),
    ],
)

op = fermionic_op.extended_jw_transformation()

a = "0000000000000111111111111111"
b = "0000000000011111111111111100"

direct = np.array(
    [
        [op.matrix_element(a, a), op.matrix_element(a, b)],
        [op.matrix_element(b, a), op.matrix_element(b, b)],
    ]
)

Hsub = fq.SubspaceHamiltonian(op, fq.Subspace([[a, b]]))

basis = np.eye(2)
matrix_free = np.column_stack([Hsub @ basis[:, i] for i in range(2)])
csr = Hsub.to_csr_linearoperator().matrix.toarray()
csr_fast = Hsub.to_csr_linearoperator_fast().matrix.toarray()

print("Direct:")
print(direct)
print("Matrix-free:")
print(matrix_free)
print("CSR:")
print(csr)
print("Fast CSR:")
print(csr_fast)

Actual output

Direct:
[[ 0.+0.j -1.+0.j]
 [-1.+0.j  0.+0.j]]

Matrix-free:
[[0. 0.]
 [0. 0.]]

CSR:
[[0. 0.]
 [0. 0.]]

Fast CSR:
[[0. 0.]
 [0. 0.]]

Expected behavior

All three SubspaceHamiltonian representations should agree with the direct matrix-element calculation.

[nonhermitian]

I believe this is because of all the specializations the type=2 side has. If one sets op.set_type(1) then the result (at least for the generated CSR matrix) is correct.

[ashsaki]

Yes, it is exactly that @nonhermitian . So, there can be spinless bitstrings that do not have clean alpha and beta separation. Major bug that I introduced in the Type 2 path.