The quantum fisher information matrix calculated by my own code is inconsistent with that by qml.metric_tensor

Hello, everyone!

I am trying to calculate the quantum fisher information matrix (QIFM) for a given circuit, and from the theory of overparameterization proposed by this paper, the QFIM should be non-full rank in case the number of parameters p>4^n (n is num of qubits). In order to verify the claim of this paper, I used qml.metric_tensor to calculate the QFIM for a given circuit with 2-qubits and 10 (depth) * 2 parameters (20 > 16=4^2), but I found that the QFIM is full rank with all eigenvalues being nonzero. This is contradictory to the overparameterization theory. The code is as follows

import pennylane as qml
import numpy as np

n_qubits = 2
cir_depth = 10
np.random.seed(1)

dev = qml.device('default.qubit', wires=n_qubits, shots=5000)

param = np.random.uniform(0.0001, 2*np.pi, cir_depth*n_qubits)

@qml.qnode(dev)

def cost_fn(weights):
    for i in range(n_qubits):
        qml.RY(np.pi/4, wires=i)
    for j in range(cir_depth):
        if j % 2 == 0:
            for i in range(n_qubits):
                qml.RX(weights[i + n_qubits*j], wires=i)
        else:
            for i in range(n_qubits):
                qml.RY(weights[i + n_qubits*j], wires=i)
        qml.CNOT(wires=[0, 1])
    return qml.expval(qml.PauliZ(0) @ qml.PauliZ(1))

metric_tensor = cost_fn.metric_tensor(param, diag_approx=False) 
eigvals, eigvecs = np.linalg.eigh(metric_tensor)

The obtained eigenvalues are as follows

    array([0.11488373, 0.12320832, 0.12535609, 0.12609588, 0.12651816,
               0.18384066, 0.18538179, 0.19106924, 0.20778353, 0.20855821,
               0.20956571, 0.23503699, 0.24635274, 0.25079696, 0.25101205,
               0.25563239, 0.25832791, 0.28002561, 0.29380703, 0.36931224])

All the eigenvalues are nonzero. I have checked the code and the claim of the over-parameterization theory about the rank of this paper serval times and I didn’t find any mistakes. So I write the code to calculate the QFIM, and for the same setting used in the case of qml.metric_tensor, the achieved eigenvalues or rank of the QFIM is significantly different from that calculated by qml.metric_tensor. So what confused me is that 'Can the QFIM calculated by qml.metric_tensor be the reasonable approximation to the true QFIM ’ or there is something wrong I haven’t found. My own code wrote with qutip is as follows
Remark: the circuit is exactly the same as the above circuit achieved by qml

import qutip as qt
import numpy as np

#tensors operators together

def genFockOp(op,position,size,levels=2,opdim=0):
    opList=[qt.qeye(levels) for x in range(size-opdim)]
    opList[position]=op
    return qt.tensor(opList)

n_qubits=2 #number qubits
depth=10 #number of layers

np.random.seed(1)

#define angles for circuit
ini_angles=np.random.uniform(0.0001, 2*np.pi, depth*n_qubits)
ini_angles = np.reshape(ini_angles, [depth, n_qubits])
n_parameters=depth*n_qubits #number of parameters of circuit

#define rotations for circuit in each layer, 1: X, 2:Y 3:Z
ini_pauli=np.zeros([depth,n_qubits], dtype=int)  # given deterministic circuit structure
for i in range(depth):
    for j in range(n_qubits):
        if i % 2 == 0:
            ini_pauli[i][j]=1
        else:
            ini_pauli[i][j]=2

#operators for circuit

levels=2#

opZ=[genFockOp(qt.sigmaz(),i,n_qubits,levels) for i in range(n_qubits)]
opX=[genFockOp(qt.sigmay(),i,n_qubits,levels) for i in range(n_qubits)]
opY=[genFockOp(qt.sigmax(),i,n_qubits,levels) for i in range(n_qubits)]
opId=genFockOp(qt.qeye(levels),0,n_qubits)

H=opZ[0]*opZ[1] #local Hamiltonian to calculate energy and gradient from

entangling_layer = qt.qip.operations.cnot(n_qubits, 0 ,1)  #CNOT

# calculate state and gradients

# list of values of gradient
gradient_list = np.zeros(n_parameters)
# save here quantum state of gradient    
grad_state_list = []


#p goes from -1 to n_parameters-1. -1 is to calculate quantum state, 0 to n_parameters rest for gradient
#calculates the quantum state U\ket{0}, as well as \partial_i U \ket{0}
#those can be used to calcute energy, gradients and the Quantum Fisher metric
for p in range(-1, n_parameters):
    counter = 0 #counter to keep back of which parameter is calculated
    initial_state=qt.tensor([qt.basis(levels,0) for i in range(n_qubits)])
    #initial layer of fixed \sqrt{H} rotations
    initial_state=qt.tensor([qt.qip.operations.ry(np.pi/4) for i in range(n_qubits)])*initial_state
    # go through depth layer
    for j in range(depth):
        rot_op = []
        #define parametrized single-qubit rotations at layer j
        for k in range(n_qubits):
            angle=ini_angles[j][k]
            if ini_pauli[j][k]==1:
                rot_op.append(qt.qip.operations.rx(angle))
            elif ini_pauli[j][k]==2:
                rot_op.append(qt.qip.operations.ry(angle))
            elif ini_pauli[j][k]==3:
                rot_op.append(qt.qip.operations.rz(angle))

            #multiply in derivative of parametrized single-qubit rotation gate at layer j for parameter of circuit p
            #this is the exact derivative
            if(counter==p):
                if(ini_pauli[j][k]==1):
                    initial_state=(-1j*opX[k]/2)*initial_state
                elif(ini_pauli[j][k]==2):
                    initial_state=(-1j*opY[k]/2)*initial_state
                elif(ini_pauli[j][k]==3):
                    initial_state=(-1j*opZ[k]/2)*initial_state

            counter+=1

        # multiply in single-qubit rotations
        initial_state=qt.tensor(rot_op)*initial_state

        # add entangling layer
        initial_state = entangling_layer * initial_state

    if p==-1: # get quantum state for cost function
        # cost function given by <\psi|H\psi}>
        circuit_state = qt.Qobj(initial_state)  # state generated by circuit
        energy = qt.expect(H, circuit_state)
        print('Energy of state', energy)

    else: # get quantum state needed for gradient
        grad_state_list.append(qt.Qobj(initial_state))  # stata with gradient applied for =-th parameter

        # gradient of circuit is given by 2*real(<\psi|H|\partial_p\psi>)
        gradient_list[p] = 2*np.real(circuit_state.overlap(H*initial_state))

# quantum fisher information metric
# calculated as \text{Re}(\braket{\partial_i \psi}{\partial_j \psi}-\braket{\partial_i \psi}{\psi}\braket{\psi}{\partial_j \psi})

# first calculate elements \braket{\psi}{\partial_j \psi}
single_qfi_elements=np.zeros(n_parameters, dtype=np.complex128)
for p in range(n_parameters):
    #print(circuit_state.overlap(grad_state_list[p]))
    single_qfi_elements[p]=circuit_state.overlap(grad_state_list[p])

# calculate the qfi matrix
qfi_matrix = np.zeros([n_parameters, n_parameters])
for p in range(n_parameters):
    for q in range(p, n_parameters):
        qfi_matrix[p][q] = np.real(grad_state_list[p].overlap(grad_state_list[q]) - single_qfi_elements[p] * single_qfi_elements[q])
        # use fact that qfi matrix is real and hermitian
        qfi_matrix[q][p]=qfi_matrix[p][q]

eigvals, eigvecs = np.linalg.eigh(qfi_matrix)

The obtained eigenvalues are as follows

array([-1.78145700e-16, -1.36750301e-16, -9.45604701e-17, -4.43704514e-17,
       -2.44116888e-18,  1.58151366e-17,  4.32981744e-17,  5.96902762e-17,
        7.05087355e-17,  1.00267602e-16,  1.60828731e-16,  3.16326547e-16,
        6.87812040e-16,  1.61709611e-01,  3.74585700e-01,  5.11669751e-01,
        6.28779511e-01,  7.28169134e-01,  1.18687949e+00,  2.67811314e+00])

The most eigenvalues calculated by my own code are closes to zero, which accords with the claim of the above-mentioned paper.

Any help will be appreciated, and someone not familiar with the package qutip can refer to this documentation or tell me in the reply, and I will give a version of numpy.

Hi @cyrie_wang! Welcome :slight_smile:

The issue here could be that PennyLane’s qml.metric_tensor function currently only supports returning the diagonal or the block-diagonal approximation to the metric tensor.

Unfortunately, this isn’t very clear in the documentation. We’re making some changes to address this:

  • In the link above, you’ll notice that the keyword argument has changed to approx="diag" or approx="block-diag", rather than True/False.

  • We’re in the process of adding support for the full metric tensor via a pull request.

But that would be my guess — if you compare qfi_matrix vs. metric_tensor , you’ll likely find that they only match on the block-diagonal.

Note: they might also differ by a factor of 2, 4, 1/2, or 1/4, depending on the convention used! The QFI/metric tensor coefficient convention is not fully settled :slight_smile:

Thanks very much for your help:grinning:, I think I probably understand what’s going on.