Add a `matrix()` for `dataset.Hamiltonian`

[minhtriet]

Feature details

When running

import pennylane as qml

h2_dataset = qml.data.load("qchem", molname="H2", bondlength=0.742, basis="STO-3G")
h2 = h2_dataset[0]
H, qubits = h2.hamiltonian, len(h2.hamiltonian.wires)

print(H.matrix())

We will have pennylane.operation.MatrixUndefinedError

Implementation

I have a working implementation that works on the above example.

Here is the Hamilton of H2

(-0.22250914236600539) [Z2]
+ (-0.22250914236600539) [Z3]
+ (-0.09963387941370971) [I0]
+ (0.17110545123720225) [Z1]
+ (0.17110545123720233) [Z0]
+ (0.12051027989546245) [Z0 Z2]
+ (0.12051027989546245) [Z1 Z3]
+ (0.16584090244119712) [Z0 Z3]
+ (0.16584090244119712) [Z1 Z2]
+ (0.16859349595532533) [Z0 Z1]
+ (0.1743207725924201) [Z2 Z3]
+ (-0.04533062254573469) [Y0 Y1 X2 X3]
+ (-0.04533062254573469) [X0 X1 Y2 Y3]
+ (0.04533062254573469) [Y0 X1 X2 Y3]
+ (0.04533062254573469) [X0 Y1 Y2 X3]

We cannot add the matrix of Z2 with Y0 Y1 X2 X3, so the main idea is for Z(2), we convert into I(0) @ Z(1) @ I(1) @ I(2), and do the same for every Hamiltonian operator in the list

Draft implementation

def hamiltonian_to_matrix(hamiltonian, wires: qml.wires.Wires):
    matrix = 0+0j
    set_wires = set(wires)
    list_wires = list(wires)
    for coeff, op in zip(hamiltonian.coeffs, hamiltonian.ops):
        missing_wires = set_wires - set(op.wires)
        padded_ops = op.obs if isinstance(op, qml.operation.Tensor) else [op]        
        for wire in list_wires:
            if wire in missing_wires:
                padded_ops.append(qml.Identity(wire))
        padded_ops.sort(key=lambda x: x.wires.tolist())
        padded_ops = [op.matrix() for op in padded_ops]        
        matrix += coeff * functools.reduce(lambda a, b: np.kron(a, b), padded_ops)
    return matrix

How important would you say this feature is?

2: Somewhat important. Needed this quarter.

Additional information

No response

[albi3ro]

Thanks for opening this issue @minhtriet . While our legacy Hamiltonian class never defined a matrix, our new analogues Sum and LinearCombination do. So you should be able to do the above on master now!

I would just be aware that computing the matrix via H.matrix() could get expensive for larger systems. In those cases, H.sparse_matrix() might be a more performant alternative if you do not need to preserve differentiability.

[DSGuala]

Additionally, qml.matrix(H) should work in this case.

On PennyLane 0.35.1, the following code prints a 16x16 matrix.

import pennylane as qml

h2_dataset = qml.data.load("qchem", molname="H2", bondlength=0.742, basis="STO-3G")
h2 = h2_dataset[0]
H, qubits = h2.hamiltonian, len(h2.hamiltonian.wires)

print(qml.matrix(H))

[albi3ro]

@minhtriet Mind if we close this issue?