How to use the pennylane.Hermitian function in PennyLane

def circuit():
            """Test QNode"""
            qml.Hadamard(wires=0)
            qml.PauliY(wires=0)
            return qml.expval(qml.Hermitian(H, 0))

def test_hermitian_diagonalizing_gates_integration(self, observable, eigvals, eigvecs, tol):
        """Tests that the diagonalizing_gates method of the Hermitian class contains contains a gate that diagonalizes the
        given observable."""
        num_wires = 2

        tensor_obs = np.kron(observable, observable)
        eigvals = np.kron(eigvals, eigvals)

        dev = qml.device('default.qubit', wires=num_wires)

        diag_gates = qml.Hermitian.diagonalizing_gates(tensor_obs, wires=list(range(num_wires)))

        assert len(diag_gates) == 1

        U = diag_gates[0].parameters[0]
        x = U @ tensor_obs @ U.conj().T
        assert np.allclose(np.diag(np.sort(eigvals)), x, atol=tol, rtol=0)

def test_hermitian_expectation(self, theta, phi, tol):
        """Test that arbitrary Hermitian expectation values are correct"""
        dev = TensorNetworkTF(wires=2)
        queue = [qml.RY(theta, wires=0), qml.RY(phi, wires=1), qml.CNOT(wires=[0, 1])]
        observables = [qml.Hermitian(A, wires=[i]) for i in range(2)]

        for i in range(len(observables)):
            observables[i].return_type = qml.operation.Expectation

        res = dev.execute(queue, observables, {})

        a = A[0, 0]
        re_b = A[0, 1].real
        d = A[1, 1]
        ev1 = ((a - d) * np.cos(theta) + 2 * re_b * np.sin(theta) * np.sin(phi) + a + d) / 2
        ev2 = ((a - d) * np.cos(theta) * np.cos(phi) + 2 * re_b * np.sin(phi) + a + d) / 2
        expected = np.array([ev1, ev2])

        assert np.allclose(res, expected, atol=tol, rtol=0)

def obs():
    """A fixture of observables to go after the queue fixture."""
    return [
        qml.expval(qml.PauliX(wires=0)),
        qml.expval(qml.Hermitian(np.identity(4), wires=[1, 2])),
    ]

def test_sample_values_hermitian(self, tol):
        """Tests if the samples of a Hermitian observable returned by sample have
        the correct values
        """
        theta = 0.543
        shots = 100000
        A = np.array([[1, 2j], [-2j, 0]])

        dev = plf.QVMDevice(device="1q-qvm", shots=shots)

        dev.apply('RX', wires=[0], par=[theta])
        dev._obs_queue = [qml.Hermitian(A, wires=[0], do_queue=False)]
        dev.pre_measure()

        s1 = dev.sample('Hermitian', [0], [A], shots)

        # s1 should only contain the eigenvalues of
        # the hermitian matrix
        eigvals = np.linalg.eigvalsh(A)
        assert np.allclose(sorted(list(set(s1))), sorted(eigvals), atol=tol, rtol=0)

        # the analytic mean is 2*sin(theta)+0.5*cos(theta)+0.5
        assert np.allclose(np.mean(s1), 2*np.sin(theta)+0.5*np.cos(theta)+0.5, atol=0.1, rtol=0)

        # the analytic variance is 0.25*(sin(theta)-4*cos(theta))^2
        assert np.allclose(np.var(s1), 0.25*(np.sin(theta)-4*np.cos(theta))**2, atol=0.1, rtol=0)

def circuit(gammas, betas, edge=None, n_layers=1):
    # apply Hadamards to get the n qubit |+> state
    for wire in range(n_wires):
        qml.Hadamard(wires=wire)
    # p instances of unitary operators
    for i in range(n_layers):
        U_C(gammas[i])
        U_B(betas[i])
    if edge is None:
        # measurement phase
        return qml.sample(comp_basis_measurement(range(n_wires)))
    # during the optimization phase we are evaluating a term
    # in the objective using expval
    return qml.expval(qml.Hermitian(pauli_z_2, wires=edge))

def layer0_off_diag_double(params):
    layer0_subcircuit(params)
    ZZ = np.kron(np.diag([1, -1]), np.diag([1, -1]))
    return expval(qml.Hermitian(ZZ, wires=[0, 1]))

for e in var_observables:
            # need to calculate d/dp
            w = e.wires

            if e.name == "Hermitian":
                # since arbitrary Hermitian observables
                # are not guaranteed to be involutory, need to take them into
                # account separately to calculate d/dp

                A = e.params[0]  # Hermitian matrix
                # if not np.allclose(A @ A, np.identity(A.shape[0])):
                new = qml.expval(qml.Hermitian(A @ A, w, do_queue=False))
            else:
                # involutory, A^2 = I
                # For involutory observables (A^2 = I) we have d/dp = 0
                new = qml.expval(qml.Hermitian(np.identity(2 ** len(w)), w, do_queue=False))

            # replace the var(A) observable with 
            self.circuit.update_node(e, new)
            new_observables.append(new)

        # calculate the analytic derivatives of the  observables
        pdA2 = self._pd_analytic(idx, args, kwargs)

        # restore the original observables, but convert their return types to expectation
        for e, new in zip(var_observables, new_observables):
            self.circuit.update_node(new, e)
            e.return_type = ObservableReturnTypes.Expectation

        # evaluate <a>
        evA = np.asarray(self.evaluate(args, kwargs))

</a>

# for each operation in the layer, get the generator and convert it to a variance
            for n, op in enumerate(curr_ops):
                gen, s = op.generator
                w = op.wires

                if gen is None:
                    raise QuantumFunctionError(
                        "Can't generate metric tensor, operation {}"
                        "has no defined generator".format(op)
                    )

                # get the observable corresponding to the generator of the current operation
                if isinstance(gen, np.ndarray):
                    # generator is a Hermitian matrix
                    variance = var(qml.Hermitian(gen, w, do_queue=False))

                    if not diag_approx:
                        Ki_matrices.append((n, expand(gen, w, self.num_wires)))

                elif issubclass(gen, Observable):
                    # generator is an existing PennyLane operation
                    variance = var(gen(w, do_queue=False))

                    if not diag_approx:
                        if issubclass(gen, qml.PauliX):
                            mat = np.array([[0, 1], [1, 0]])
                        elif issubclass(gen, qml.PauliY):
                            mat = np.array([[0, -1j], [1j, 0]])
                        elif issubclass(gen, qml.PauliZ):
                            mat = np.array([[1, 0], [0, -1]])

def prepare_and_sample(weights):

    # Variational circuit generating a guess for the solution vector |x>
    variational_block(weights)

    # We assume that the system is measured in the computational basis.
    # If we label each basis state with a decimal integer j = 0, 1, ... 2 ** n_qubits - 1,
    # this is equivalent to a measurement of the following diagonal observable.
    basis_obs = qml.Hermitian(np.diag(range(2 ** n_qubits)), wires=range(n_qubits))

    return qml.sample(basis_obs)