Numerically computing fidelities close to 1

TL;DR

For two states \(\ket{a}\) and \(\ket{b}\), let \(c = \braket{a|b}/|\braket{a|b}|\) and \(e = |c\ket{a} - \ket{b}|^2\) (the norm squared of the difference vector between \(c\ket{a}\) and \(\ket{b}\)). Then the fidelity \(\mathcal{F} = |\braket{a|b}|^2\) can be computed with good numerical precision as \[1 - \mathcal{F} = e - e^2/4\]

Note

I doubt I’m the first person to come up with this trick, but I wasn’t able to find it anywhere on the internet when I needed it. Let me know if you know a reference for it!

Suppose we want to examine the fidelity \(\mathcal{F} = |\braket{a|b}|^2\) of two quantum states \(\ket{a}\) and \(\ket{b}\), stored as vectors of complex numbers. If \(\ket{a}\) and \(\ket{b}\) are very close, it is more convenient to study the “error” \(\epsilon = 1-\mathcal{F}\) rather than the fidelity itself. However, computing this value in the obvious way runs into problems with floating-point precision. Here’s a straightforward function computing \(1-\mathcal{F}\):

def error_naive(a, b):
    return 1 - np.abs(np.vdot(a, b))**2

Let’s try plotting the measured fidelity vs. actual fidelity, for a range of fidelities very close to 1:

Code

import numpy as np
from matplotlib import pyplot as plt

def generate_random_complex_unit_vector(dimension):
    rtn = np.zeros((dimension,), dtype=np.complex64)
    rtn += np.random.standard_normal(dimension)
    rtn += 1j * np.random.standard_normal(dimension)
    rtn /= np.linalg.norm(rtn)
    return rtn

def generate_random_orthogonal_vectors(dimension):
    a = generate_random_complex_unit_vector(dimension)
    b = generate_random_complex_unit_vector(dimension)

    # gram-schmidt, to make them orthogonal
    b -= np.vdot(a, b) * a
    b /= np.linalg.norm(b)

    return a, b

def plot_error_measure(dimension, x_min, x_max, n_pts, error_funcs):
    a, b = generate_random_orthogonal_vectors(dimension)
    
    pts = np.geomspace(x_min, x_max, n_pts)

    for error_func in error_funcs:
        results = []
        for target_err in pts:
            c = np.sqrt(target_err)*b + np.sqrt(1-target_err)*a
            results.append(error_func(a, c))

        plt.loglog(pts, results, label=error_func.__name__)
    
    plt.loglog(pts, pts, color='k', linestyle=':', label=r'$1 - \mathcal{F}$')

    plt.xlabel(r'$1-\mathcal{F}$')
    plt.ylabel('Result')

    plt.title('Computed vs. actual fidelity')

    plt.legend()

plot_error_measure(
    dimension=2**20, 
    x_min=1E-14, 
    x_max=1, 
    n_pts=21, 
    error_funcs=[error_naive]
)

The reason for the discrepancy is that in the inner product, we are summing values \(a_i^\dagger b_i\) that have magnitude \(\mathcal{O}(1/d)\), but the “important part” of each term (that prevents the sum from being exactly 1) is of order \(\mathcal{O}(\epsilon/d)\). If \(\epsilon\) is very small, this “content” will be lost to floating point errors. Instead, we’d like to figure out a way to sum that content we care about directly.

Consider the following fact: \[ |\ket{a}-\ket{b}|^2 = 2 - 2\mathrm{Re}\left[ \braket{a | b} \right] \] So, the square of the norm of the difference gives the ~vibes~ of \(\sqrt{\epsilon} = |\braket{a|b}|\), except it cares about the relative phase between \(\ket{a}\) and \(\ket{b}\). The upside is that it should be numerically stable to compute, because the elements of the difference vector will be of order \(\mathcal{O}(\epsilon/d)\). The key idea is that we can use it to measure \(\epsilon\) if we ensure there is no relative phase between \(\ket{a}\) and \(\ket{b}\).

Let \(c = \braket{a|b}/|\braket{a|b}|\). Then by construction \(\braket{c a | b} = |\braket{a | b}| = \mathrm{Re}\left[ \braket{c a | b} \right]\). Combining this with the above, we may define \(e = |c \ket{a}-\ket{b}|^2\) and then: \[ |\braket{a | b}| = 1 - e/2 \] and thus \[ 1-\mathcal{F} = 1 - |\braket{a|b}|^2 = e - e^2/4 \] This is an exact expression (we didn’t use any approximations to get it), but it also will not suffer from numerical instability for values close to zero.

Pretty cool! Let’s try it:

def error_precise(a, b):
    inner_prod = np.vdot(a, b)

    # to avoid division by zero
    if inner_prod == 0:
        return 1.0
        
    c = inner_prod / np.abs(inner_prod)
    diff = c*a - b
    e = np.real(np.vdot(diff, diff))  # imag part is zero
    return e - e**2 / 4

Code

plot_error_measure(
    dimension=2**20, 
    x_min=1E-14, 
    x_max=1, 
    n_pts=21, 
    error_funcs=[error_naive, error_precise]
)

It works super well! Yay!