Fidelity and random matrices

docs/code.rst

@@ -14,6 +14,8 @@ Python interface
 
 * The :mod:`thewalrus.symplectic` submodule provides access to a convenient set of symplectic matrices and utility functions to manipulate them
 
+* The :mod:`thewalrus.random` submodule provides access to random unitary, symplectic and covariance matrices
+
 * The :mod:`thewalrus.fock_gradients` submodule provides access to the Fock representation of certain continuous-variable gates and their gradients
 
 * The :mod:`thewalrus.reference` submodule provides access to pure-Python reference implementations of the hafnian, loop hafnian, and torontonian

docs/code/random.rst

@@ -0,0 +1,2 @@
+.. automodule:: thewalrus.random
+    :members:

docs/gbs_sampling.rst

@@ -80,9 +80,9 @@ In the previous section it was mentioned that states whose covariance matrix sat
 
    To generate photon number or threshold samples from a classical gaussian state specified by a quadrature covariance matrix use :func:`thewalrus.samples.hafnian_sample_classical_state` or :func:`thewalrus.samples.torontonian_sample_classical_state`.
 
-Note that one can use this observation to sample from a probability distribution that is proportional to the permanent of positive semidefinite matrix, for detail of how this is done cf. Ref. :cite:`jahangiri2020point`.
+Note that one can use this observation to sample from a probability distribution that is proportional to the permanent of a positive semidefinite matrix, for details on how this is done cf. Ref. :cite:`jahangiri2020point`.
 
 .. tip::
 
-   To generate photon number samples from thermal state parametrized by a positive semidefinite *real* matrix use the module :func:`thewalrus.csamples`.
+   To generate photon number samples from a thermal state parametrized by a positive semidefinite *real* matrix use the module :func:`thewalrus.csamples`.
 

docs/index.rst

@@ -165,6 +165,7 @@ The Walrus library is **free** and **open source**, released under the Apache Li
    code/samples
    code/csamples
    code/symplectic
+   code/random
    code/fock_gradients
    code/reference
    code/libwalrus

thewalrus/quantum.py

@@ -53,6 +53,7 @@
     probabilities
     update_probabilities_with_loss
     update_probabilities_with_noise
+    fidelity
 
 Details
 ^^^^^^^
@@ -120,6 +121,7 @@
 from scipy.optimize import root_scalar
 from scipy.special import factorial as fac
 from scipy.stats import nbinom
+from scipy.linalg import sqrtm
 from numba import jit
 
 from thewalrus.symplectic import expand, sympmat, is_symplectic
@@ -1232,3 +1234,55 @@ def update_probabilities_with_noise(probs_noise, probs):
         probs_masked = _update_1d(probs_masked, one_d, cutoff)
         probs = np.transpose(probs_masked, axes=perm)
     return probs
+
+
+def fidelity(mu1, mu2, cov1, cov2, hbar=2, rtol=1e-05, atol=1e-08):
+    """Calculates the fidelity between two Gaussian quantum states.
+    Note that if the covariance matrices correspond to pure states this
+    function reduces to the modulus square of the overlap of their state vectors.
+    For the derivation see:
+
+    `'Quantum Fidelity for Arbitrary Gaussian States', Banchi et al. <10.1103/PhysRevLett.115.260501>`_
+
+    Args:
+        mu1 (array): vector of means of the first state
+        mu2 (array): vector of means of the second state
+        cov1 (array): covariance matrix of the first state
+        cov2 (array): covariance matrix of the second state
+        hbar (float): value of hbar in the uncertainty relation
+        rtol (float): the relative tolerance parameter used in `np.allclose`
+        atol (float): the absolute tolerance parameter used in `np.allclose`
+
+    Returns:
+        (float): value of the fidelity between the two states
+    """
+
+    n0, n1 = cov1.shape
+    m0, m1 = cov2.shape
+    (l0,) = mu1.shape
+    (l1,) = mu1.shape
+    if not n0 == n1 == m0 == m1 == l0 == l1:
+        raise ValueError("The inputs have incompatible shapes")
+
+    v1 = cov1 / hbar
+    v2 = cov2 / hbar
+    deltar = (mu1 - mu2) / np.sqrt(hbar / 2)
+    n0, n1 = cov1.shape
+    n = n0 // 2
+    W = sympmat(n)
+
+    si12 = np.linalg.inv(v1 + v2)
+    vaux = W.T @ si12 @ (0.25 * W + v2 @ W @ v1)
+    p1 = vaux @ W
+    p1 = p1 @ p1
+    p1 = np.identity(2 * n) + 0.25 * np.linalg.inv(p1)
+    if np.allclose(p1, 0, rtol=rtol, atol=atol):
+        p1 = np.zeros_like(p1)
+    else:
+        p1 = sqrtm(p1)
+    p1 = 2 * (p1 + np.identity(2 * n))
+    p1 = p1 @ vaux
+    f = np.sqrt(np.linalg.det(si12) * np.linalg.det(p1)) * np.exp(
+        -0.25 * deltar @ si12 @ deltar
+    )
+    return f

thewalrus/random.py

@@ -0,0 +1,118 @@
+# Copyright 2019 Xanadu Quantum Technologies Inc.
+
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+#     http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+"""
+Random matrices
+===============
+
+.. currentmodule:: thewalrus.random
+
+This submodule provides access to utility functions to generate random unitary, symplectic
+and covariance matrices.
+
+"""
+import numpy as np
+import scipy as sp
+
+# ------------------------------------------------------------------------
+# Random numbers and matrices                                           |
+# ------------------------------------------------------------------------
+
+def randnc(*arg):
+    """Normally distributed array of random complex numbers."""
+    return np.random.randn(*arg) + 1j * np.random.randn(*arg)
+
+
+def random_covariance(N, hbar=2, pure=False, block_diag=False):
+    r"""Random covariance matrix.
+
+    Args:
+        N (int): number of modes
+        hbar (float): the value of :math:`\hbar` to use in the definition
+            of the quadrature operators :math:`x` and :math:`p`
+        pure (bool): If True, a random covariance matrix corresponding
+            to a pure state is returned.
+        block_diag (bool): If True, uses passive Gaussian transformations that are orthogonal
+            instead of unitary. This implies that the positions :math:`x` do not mix with
+            the momenta :math:`p` and thus the covariance matrix is block diagonal.
+
+    Returns:
+        array: random :math:`2N\times 2N` covariance matrix
+    """
+    S = random_symplectic(N, block_diag=block_diag)
+
+    if pure:
+        return (hbar / 2) * S @ S.T
+
+    nbar = 2 * np.abs(np.random.random(N)) + 1
+    Vth = (hbar / 2) * np.diag(np.concatenate([nbar, nbar]))
+
+    return S @ Vth @ S.T
+
+
+def random_symplectic(N, passive=False, block_diag=False, scale=1.0):
+    r"""Random symplectic matrix representing a Gaussian transformation.
+    The squeezing parameters :math:`r` for active transformations are randomly
+    sampled from the standard normal distribution, while passive transformations
+    are randomly sampled from the Haar measure. Note that for the Symplectic
+    group there is no notion of Haar measure since this is group is not compact.
+
+    Args:
+        N (int): number of modes
+        passive (bool): If True, returns a passive Gaussian transformation (i.e.,
+            one that preserves photon number). If False (default), returns an active
+            transformation.
+        block_diag (bool): If True, uses passive Gaussian transformations that are orthogonal
+            instead of unitary. This implies that the positions :math:`q` do not mix with
+            the momenta :math:`p` and thus the symplectic operator is block diagonal
+        scale (float): Sets the scale of the random values used as squeezing parameters.
+            They will range from 0 to :math:`\sqrt{2}\texttt{scale}`
+
+    Returns:
+        array: random :math:`2N\times 2N` symplectic matrix
+    """
+    U = random_interferometer(N, real=block_diag)
+    O = np.block([[U.real, -U.imag], [U.imag, U.real]])
+
+    if passive:
+        return O
+
+    U = random_interferometer(N, real=block_diag)
+    P = np.block([[U.real, -U.imag], [U.imag, U.real]])
+
+    r = scale * np.abs(randnc(N))
+    Sq = np.diag(np.concatenate([np.exp(-r), np.exp(r)]))
+
+    return O @ Sq @ P
+
+
+def random_interferometer(N, real=False):
+    r"""Random unitary matrix representing an interferometer.
+    For more details, see :cite:`mezzadri2006`.
+
+    Args:
+        N (int): number of modes
+        real (bool): return a random real orthogonal matrix
+
+    Returns:
+        array: random :math:`N\times N` unitary distributed with the Haar measure
+    """
+    if real:
+        z = np.random.randn(N, N)
+    else:
+        z = randnc(N, N) / np.sqrt(2.0)
+    q, r = sp.linalg.qr(z)
+    d = np.diagonal(r)
+    ph = d / np.abs(d)
+    U = np.multiply(q, ph, q)
+    return U

thewalrus/samples.py

@@ -483,8 +483,8 @@ def torontonian_sample_graph(A, n_mean, samples=1, max_photons=30, parallel=Fals
 def hafnian_sample_classical_state(
     cov, samples, mean=None, hbar=2, atol=1e-08, cutoff=None
 ):  # add cutoff for consistency pylint: disable=unused-argument
-    r"""Returns samples from a Gaussian state that has a positive :math:`P` function.
 
+    r"""Returns samples from a Gaussian state that has a positive :math:`P` function.
     Args:
         cov(array): a :math:`2N\times 2N` ``np.float64`` covariance matrix
             representing an :math:`N` mode quantum state. This can be obtained

thewalrus/tests/test_quantum.py

@@ -18,11 +18,12 @@
 import pytest
 
 import numpy as np
-from scipy.linalg import qr
 from scipy.stats import poisson
 
 from thewalrus.symplectic import rotation, squeezing, interferometer, two_mode_squeezing, beam_splitter, loss
 
+from thewalrus.random import random_covariance, random_interferometer
+
 from thewalrus.quantum import (
     reduced_gaussian,
     Xmat,
@@ -53,36 +54,14 @@
     update_probabilities_with_loss,
     update_probabilities_with_noise,
     loss_mat,
+    fidelity,
 )
 
 
 # make tests deterministic
 np.random.seed(137)
 
 
-def random_interferometer(N, real=False):
-    r"""Random unitary matrix representing an interferometer.
-
-    For more details, see arXiv:math-ph/0609050
-
-    Args:
-        N (int): number of modes
-        real (bool): return a random real orthogonal matrix
-
-    Returns:
-        array: random :math:`N\times N` unitary distributed with the Haar measure
-    """
-    if real:
-        z = np.random.randn(N, N)
-    else:
-        z = (np.random.randn(N, N) + 1j * np.random.randn(N, N)) / np.sqrt(2.0)
-    q, r = qr(z)
-    d = np.diag(r)
-    ph = d / np.abs(d)
-    U = np.multiply(q, ph, q)
-    return U
-
-
 @pytest.mark.parametrize("n", [0, 1, 2])
 def test_reduced_gaussian(n):
     """test that reduced gaussian returns the correct result"""
@@ -1042,3 +1021,84 @@ def test_update_with_noise_coherent_value_error():
         match="The list of probability distributions probs_noise and the tensor of probabilities probs have incompatible dimensions.",
     ):
         update_probabilities_with_noise(noise_dists, probs)
+
+
+@pytest.mark.parametrize("hbar", [1 / 2, 1, 2, 1.6])
+@pytest.mark.parametrize("num_modes", np.arange(5, 10))
+@pytest.mark.parametrize("pure", [True, False])
+@pytest.mark.parametrize("block_diag", [True, False])
+def test_fidelity_with_self(num_modes, hbar, pure, block_diag):
+    """Test that the fidelity of two identical quantum states is 1"""
+    cov = random_covariance(num_modes, hbar=hbar, pure=pure, block_diag=block_diag)
+    means = np.random.rand(2 * num_modes)
+    assert np.allclose(fidelity(means, means, cov, cov, hbar=hbar), 1, atol=1e-4)
+
+
+@pytest.mark.parametrize("hbar", [1 / 2, 1, 2, 1.6])
+@pytest.mark.parametrize("num_modes", np.arange(5, 10))
+@pytest.mark.parametrize("pure", [True, False])
+@pytest.mark.parametrize("block_diag", [True, False])
+def test_fidelity_is_symmetric(num_modes, hbar, pure, block_diag):
+    """Test that the fidelity is symmetric and between 0 and 1"""
+    cov1 = random_covariance(num_modes, hbar=hbar, pure=pure, block_diag=block_diag)
+    means1 = np.sqrt(2 * hbar) * np.random.rand(2 * num_modes)
+    cov2 = random_covariance(num_modes, hbar=hbar, pure=pure, block_diag=block_diag)
+    means2 = np.sqrt(2 * hbar) * np.random.rand(2 * num_modes)
+    f12 = fidelity(means1, means2, cov1, cov2, hbar=hbar)
+    f21 = fidelity(means2, means1, cov2, cov1, hbar=hbar)
+    assert np.allclose(f12, f21)
+    assert 0 <= np.real_if_close(f12) < 1.0
+
+
+@pytest.mark.parametrize("num_modes", np.arange(5, 10))
+@pytest.mark.parametrize("hbar", [0.5, 1, 2, 1.6])
+def test_fidelity_coherent_state(num_modes, hbar):
+    """Test the fidelity of two multimode coherent states"""
+    beta1 = np.random.rand(num_modes) + 1j * np.random.rand(num_modes)
+    beta2 = np.random.rand(num_modes) + 1j * np.random.rand(num_modes)
+    means1 = Means(np.concatenate([beta1, beta1.conj()]), hbar=hbar)
+    means2 = Means(np.concatenate([beta2, beta2.conj()]), hbar=hbar)
+    cov1 = hbar * np.identity(2 * num_modes) / 2
+    cov2 = hbar * np.identity(2 * num_modes) / 2
+    fid = fidelity(means1, means2, cov1, cov2, hbar=hbar)
+    expected = np.exp(-np.linalg.norm(beta1 - beta2) ** 2)
+    assert np.allclose(expected, fid)
+
+
+@pytest.mark.parametrize("hbar", [0.5, 1, 2, 1.6])
+@pytest.mark.parametrize("r", [-2, 0, 2])
+@pytest.mark.parametrize("alpha", np.random.rand(10) + 1j * np.random.rand(10))
+def test_fidelity_vac_to_displaced_squeezed(r, alpha, hbar):
+    """Calculates the fidelity between a coherent squeezed state and vacuum"""
+    cov1 = np.diag([np.exp(2 * r), np.exp(-2 * r)]) * hbar / 2
+    means1 = Means(np.array([alpha, np.conj(alpha)]), hbar=hbar)
+    means2 = np.zeros([2])
+    cov2 = np.identity(2) * hbar / 2
+    expected = (
+        np.exp(-np.abs(alpha) ** 2)
+        * np.abs(np.exp(np.tanh(r) * np.conj(alpha) ** 2))
+        / np.cosh(r)
+    )
+    assert np.allclose(expected, fidelity(means1, means2, cov1, cov2, hbar=hbar))
+
+
+@pytest.mark.parametrize("hbar", [0.5, 1, 2, 1.6])
+@pytest.mark.parametrize("r1", np.random.rand(3))
+@pytest.mark.parametrize("r2", np.random.rand(3))
+def test_fidelity_squeezed_vacuum(r1, r2, hbar):
+    """Tests fidelity between two squeezed states"""
+    cov1 = np.diag([np.exp(2 * r1), np.exp(-2 * r1)]) * hbar / 2
+    cov2 = np.diag([np.exp(2 * r2), np.exp(-2 * r2)]) * hbar / 2
+    mu = np.zeros([2])
+    assert np.allclose(1 / np.cosh(r1 - r2), fidelity(mu, mu, cov1, cov2, hbar=hbar))
+
+
+def test_fidelity_wrong_shape():
+    """Tests the correct error is raised"""
+    cov1 = np.identity(2)
+    cov2 = np.identity(4)
+    mu = np.zeros(2)
+    with pytest.raises(
+        ValueError, match="The inputs have incompatible shapes"
+    ):
+        fidelity(mu, mu, cov1, cov2)