Optimizing state overlap / multi-qubit qml.expval.Hermitian()

And just for completeness, here is the full optimization:

import pennylane as qml
from pennylane import numpy as np

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

@qml.qnode(dev)
def circuit(x, target_observable=None):
    qml.RX(x[0], wires=0)
    qml.RY(x[1], wires=0)
    qml.RZ(x[2], wires=0)
    qml.CNOT(wires=[0, 1])
    return qml.expval.Hermitian(target_observable, wires=[0, 1])

target_state = 1/np.sqrt(2) * np.array([1,0,0,1])
target_herm_op = np.outer(target_state.conj(), target_state)

weights = np.array([0.5, 0.1, 0.2])

def cost(weights):
    return np.abs(circuit(weights, target_observable=target_herm_op)-1)**2

opt = qml.AdamOptimizer(stepsize=0.4)

for i in range(2000):
    weights = opt.step(cost, weights)

    if (i+1) % 20 == 0:
        print('Cost after step {:5d}: {: .7f}'.format(i+1, cost(weights)))

fidelity = circuit(weights, target_observable=target_herm_op)

print('Optimized rotation angles: {}'.format(weights))
print('Final fidelity: {}'.format(fidelity))

# this is a hidden/undocumented attribute of the default.qubit
# plugin, that is not part of the API.
print('Output state: {}'.format(dev._state))

which produces the output

Optimized rotation angles: [ 1.44393011  1.57000406  1.44393011]
Final fidelity: 0.9999999974878844
Output state: [ 0.70714217-0.00027786j,  0.00000000+0.j  0.00000000+0.j  0.70707128-0.00027786j]