Fidelity calculation for density matrices improvement

[alvinquantum]

Cirq/cirq-core/cirq/qis/measures.py

Lines 239 to 242 in 782c14e

	state1_sqrt = _sqrt_positive_semidefinite_matrix(state1)
	eigs = scipy.linalg.eigvalsh(state1_sqrt @ state2 @ state1_sqrt)
	trace = np.sum(np.sqrt(np.abs(eigs)))
	return trace ** 2

The current implementation uses the direct definition of fidelity, which works well for small dimensions. From my experience, I was getting values significantly greater than 1. For the fidelity between rho1 and rho2, I found that squaring the sum of the singular values of the product of square_root(rho1) and square_root(rho2) gives a more accurate calculation. This equation is equal to fidelity. The idea is not mine and came from qutip.

This is my version of the implementation that I use.

import scipy
def get_fidelity(rho1, rho2):
    rho1_sqrt=scipy.linalg.sqrtm(rho1)
    rho2_sqrt=scipy.linalg.sqrtm(rho2)
    return scipy.linalg.svdvals(rho1_sqrt @ rho2_sqrt).sum()**2

[tanujkhattar]

This looks like a reasonable request to me.

@viathor What do you think?

For completeness, the modified formula can be derived as follows:

• F(rho1, rho2) = tr(|rho1_sqrt @ rho2_sqrt|) ** 2, where |A| = sqrt(A.T @ A)
  • Ref: https://en.wikipedia.org/wiki/Fidelity_of_quantum_states
• tr(|A|) = scipy.linalg.svdvals(A).sum()
  • Ref: https://math.stackexchange.com/questions/1840742/a-confusion-about-the-definition-of-the-trace-norm

[viathor]

I agree the request is reasonable.

Ideally, the PR implementing this would include a unit test demonstrating improved numerical stability. It should fail for the current implementation and succeed for the new one. @alvinquantum, any of the cases you encountered where fidelity exceeded 1 in cirq and stayed <=1 in your code can easily become such a unit test. Also, it'd be nice to have some evaluation of the effect on performance, e.g. ipython's %timeit output posted here would work.

@tanujkhattar Your derivation LGTM with a nit that |A| = sqrt(A.T.conj() @ A).

[alvinquantum]

I think this is a working example. Here cirq.fidelity gives around 1.09 while get_fidelity gives 1 with a 10^-5 tolerance. The code starts with a 10 qubit state, traces out the 10th system, and then calculates the fidelity of the reduced state with itself. The get_fidelity can also be improved since it tends to give values greater than 1, but from my experience those values have been greater by around 10^-5.

from cirq.contrib.qasm_import import circuit_from_qasm
import scipy
import cirq
import numpy as np
from cirq import partial_trace, DensityMatrixSimulator

def get_fidelity(rho1, rho2):
    '''Fidelity through singular values of sqrt(rho1)*sqrt(rho2)'''
    rho1_sqrt=scipy.linalg.sqrtm(rho1)
    rho2_sqrt=scipy.linalg.sqrtm(rho2)
    return scipy.linalg.svdvals(rho1_sqrt @ rho2_sqrt).sum()**2

if __name__ == "__main__":
    circ = circuit_from_qasm(
                    '''OPENQASM 2.0;
                    include "qelib1.inc";
                    qreg q[10];
                    rx(1.0349347) q[0];
                    ry(1.1344214) q[0];
                    rz(0.49872878) q[0];
                    rx(1.2483498) q[1];
                    ry(1.5452867) q[1];
                    rz(1.1293869) q[1];
                    rx(1.3588455) q[2];
                    ry(0.49602024) q[2];
                    rz(0.21943984) q[2];
                    rx(1.9583295) q[3];
                    ry(1.8588403) q[3];
                    rz(1.2323985) q[3];
                    rx(0.027308104) q[4];
                    ry(1.0267784) q[4];
                    rz(1.0619419) q[4];
                    rx(0.89107072) q[5];
                    ry(0.86430418) q[5];
                    rz(0.98182118) q[5];
                    rx(0.97120279) q[6];
                    ry(1.5382162) q[6];
                    rz(1.9003976) q[6];
                    rx(0.34264218) q[7];
                    ry(0.21962922) q[7];
                    rz(1.6263407) q[7];
                    rx(0.62256334) q[8];
                    ry(0.47124841) q[8];
                    rz(0.61290254) q[8];
                    rx(0.0017523425) q[9];
                    ry(1.7104261) q[9];
                    rz(0.13981947) q[9];
                    ''')

    init_qubits=10
    final_qubits=init_qubits-1
    keep_idxs=list(range(final_qubits))
    simulator=DensityMatrixSimulator()
    rho=simulator.simulate(circ).final_density_matrix
    rho_reshaped=np.reshape(rho, (2,2)*init_qubits)
    rho_reduced=np.reshape(partial_trace(rho_reshaped, keep_idxs),(2**final_qubits, 2**final_qubits))
    cirq_fidelity=cirq.fidelity(rho_reduced, rho_reduced, validate=False, qid_shape=(2,)*final_qubits)
    print(f"cirq.fidelity: {cirq_fidelity}") #approx. 1.0892252935292532
    print(f"get_fidelity: {get_fidelity(rho_reduced, rho_reduced)}") #approx 1.0000061793639072

[viathor]

Thanks, @alvinquantum! This is helpful.

Note: this is good basis for a unit test, but should use cirq gates rather than converting from qasm (because ideally the unit test shouldn't be sensitive to bugs in our qasm converter!).

[anonymousr007]

Hi Google Quantum Team, Is this issue still open? If, yes please tell me more about it.

[tanujkhattar]

@anonymousr007 Yes, the issue is still open.

The goal is to replace the current implementation of calculating fidelity of two density matrices with the improved implementation that has lower numerical error.

Current implementation:

Cirq/cirq-core/cirq/qis/measures.py

Lines 239 to 242 in 782c14e

	state1_sqrt = _sqrt_positive_semidefinite_matrix(state1)
	eigs = scipy.linalg.eigvalsh(state1_sqrt @ state2 @ state1_sqrt)
	trace = np.sum(np.sqrt(np.abs(eigs)))
	return trace ** 2

New implementation:

def get_fidelity(rho1, rho2):
    rho1_sqrt=scipy.linalg.sqrtm(rho1)
    rho2_sqrt=scipy.linalg.sqrtm(rho2)
    return scipy.linalg.svdvals(rho1_sqrt @ rho2_sqrt).sum()**2

The PR implementing the change should also add a unit test demonstrating improved numerical stability of the new implementation, s.t. the current implementation fails but the new one passes.

[anonymousr007]

Please, assign this issue to me. I want to work on this.