Source code for thewalrus.fock_gradients

# 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.
"""
Fock representations
====================

**Module name:** :mod:`thewalrus.fock_gradients`

.. currentmodule:: thewalrus.fock_gradients

This module contains the Fock representation of the standard Gaussian gates as well as their gradients.

Summary
-------

.. autosummary::
	displacement
	squeezing
	beamsplitter
	two_mode_squeezing
	mzgate
	grad_displacement
	grad_squeezing
	grad_beamsplitter
	grad_two_mode_squeezing
	grad_mzgate

Code details
------------
"""
import numpy as np
from numba import jit

SQRT = np.sqrt(np.arange(1e6))


[docs]@jit(nopython=True) def displacement(r, phi, cutoff, dtype=np.complex128): # pragma: no cover r"""Calculates the matrix elements of the displacement gate using a recurrence relation. Uses the log of the matrix elements to avoid numerical issues and then takes the exponential. Args: r (float): displacement magnitude phi (float): displacement angle cutoff (int): Fock ladder cutoff dtype (data type): Specifies the data type used for the calculation Returns: array[complex]: matrix representing the displacement operation. """ D = np.zeros((cutoff, cutoff), dtype=dtype) rng = np.arange(cutoff) rng[0] = 1 log_k_fac = np.cumsum(np.log(rng)) for n_minus_m in range(cutoff): m_max = cutoff - n_minus_m logL = np.log(_laguerre(r**2.0, m_max, n_minus_m)) for m in range(m_max): n = n_minus_m + m D[n, m] = np.exp( 0.5 * (log_k_fac[m] - log_k_fac[n]) + n_minus_m * np.log(r) - (r**2.0) / 2.0 + 1j * phi * n_minus_m + logL[m] ) D[m, n] = (-1.0) ** (n_minus_m) * np.conj(D[n, m]) return D
@jit(nopython=True, cache=True) def _laguerre(x, N, alpha, dtype=np.complex128): # pragma: no cover r"""Returns the N first generalized Laguerre polynomials evaluated at x. Args: x (float): point at which to evaluate the polynomials N (int): maximum Laguerre polynomial to calculate alpha (float): continuous parameter for the generalized Laguerre polynomials """ L = np.zeros(N, dtype=dtype) L[0] = 1.0 if N > 1: for m in range(0, N - 1): L[m + 1] = ((2 * m + 1 + alpha - x) * L[m] - (m + alpha) * L[m - 1]) / (m + 1) return L
[docs]@jit(nopython=True) def grad_displacement(T, r, phi): # pragma: no cover r"""Calculates the gradients of the displacement gate with respect to the displacement magnitude and angle. Args: T (array[complex]): array representing the gate r (float): displacement magnitude phi (float): displacement angle Returns: tuple[array[complex], array[complex]]: The gradient of the displacement gate with respect to r and phi """ cutoff = T.shape[0] dtype = T.dtype ei = np.exp(1j * phi) eic = np.exp(-1j * phi) alpha = r * ei alphac = r * eic sqrt = np.sqrt(np.arange(cutoff, dtype=dtype)) grad_r = np.zeros((cutoff, cutoff), dtype=dtype) grad_phi = np.zeros((cutoff, cutoff), dtype=dtype) for m in range(cutoff): for n in range(cutoff): grad_r[m, n] = -r * T[m, n] + sqrt[m] * ei * T[m - 1, n] - sqrt[n] * eic * T[m, n - 1] grad_phi[m, n] = ( sqrt[m] * 1j * alpha * T[m - 1, n] + sqrt[n] * 1j * alphac * T[m, n - 1] ) return grad_r, grad_phi
[docs]@jit(nopython=True) def squeezing(r, theta, cutoff, dtype=np.complex128): # pragma: no cover r"""Calculates the matrix elements of the squeezing gate using a recurrence relation. Args: r (float): squeezing magnitude theta (float): squeezing angle cutoff (int): Fock ladder cutoff dtype (data type): Specifies the data type used for the calculation Returns: array[complex]: matrix representing the squeezing gate. """ S = np.zeros((cutoff, cutoff), dtype=dtype) sqrt = np.sqrt(np.arange(cutoff, dtype=dtype)) eitheta_tanhr = np.exp(1j * theta) * np.tanh(r) sechr = 1.0 / np.cosh(r) R = np.array( [ [-eitheta_tanhr, sechr], [sechr, np.conj(eitheta_tanhr)], ] ) S[0, 0] = np.sqrt(sechr) for m in range(2, cutoff, 2): S[m, 0] = sqrt[m - 1] / sqrt[m] * R[0, 0] * S[m - 2, 0] for m in range(0, cutoff): for n in range(1, cutoff): if (m + n) % 2 == 0: S[m, n] = ( sqrt[n - 1] / sqrt[n] * R[1, 1] * S[m, n - 2] + sqrt[m] / sqrt[n] * R[0, 1] * S[m - 1, n - 1] ) return S
[docs]@jit(nopython=True) def grad_squeezing(T, r, phi): # pragma: no cover r"""Calculates the gradients of the squeezing gate with respect to the squeezing magnitude and angle Args: T (array[complex]): array representing the gate r (float): squeezing magnitude phi (float): squeezing angle Returns: tuple[array[complex], array[complex]]: The gradient of the squeezing gate with respect to the r and phi """ cutoff = T.shape[0] dtype = T.dtype sqrt = np.sqrt(np.arange(cutoff, dtype=dtype)) grad_r = np.zeros((cutoff, cutoff), dtype=dtype) grad_phi = np.zeros((cutoff, cutoff), dtype=dtype) sechr = 1.0 / np.cosh(r) tanhr = np.tanh(r) eiphi = np.exp(1j * phi) eiphiconj = np.exp(-1j * phi) for m in range(cutoff): for n in range(cutoff): grad_r[m, n] = ( -0.5 * tanhr * T[m, n] - sechr * tanhr * sqrt[m] * sqrt[n] * T[m - 1, n - 1] - 0.5 * eiphi * sechr**2 * sqrt[m] * sqrt[m - 1] * T[m - 2, n] + 0.5 * eiphiconj * sechr**2 * sqrt[n] * sqrt[n - 1] * T[m, n - 2] ) grad_phi[m, n] = ( -0.5j * eiphi * tanhr * sqrt[m] * sqrt[m - 1] * T[m - 2, n] - 0.5j * eiphiconj * tanhr * sqrt[n] * sqrt[n - 1] * T[m, n - 2] ) return grad_r, grad_phi
[docs]@jit(nopython=True) def two_mode_squeezing(r, theta, cutoff, dtype=np.complex128): # pragma: no cover """Calculates the matrix elements of the two-mode squeezing gate recursively. Args: r (float): squeezing magnitude theta (float): squeezing angle cutoff (int): Fock ladder cutoff dtype (data type): Specifies the data type used for the calculation Returns: array[float]: The Fock representation of the gate """ sqrt = np.sqrt(np.arange(cutoff, dtype=dtype)) Z = np.zeros((cutoff, cutoff, cutoff, cutoff), dtype=dtype) sc = 1.0 / np.cosh(r) eiptr = np.exp(-1j * theta) * np.tanh(r) R = -np.array( [ [0, -np.conj(eiptr), -sc, 0], [-np.conj(eiptr), 0, 0, -sc], [-sc, 0, 0, eiptr], [0, -sc, eiptr, 0], ] ) Z[0, 0, 0, 0] = sc # rank 2 for n in range(1, cutoff): Z[n, n, 0, 0] = R[0, 1] * Z[n - 1, n - 1, 0, 0] # rank 3 for m in range(cutoff): for n in range(m): p = m - n if 0 < p < cutoff: Z[m, n, p, 0] = R[0, 2] * sqrt[m] / sqrt[p] * Z[m - 1, n, p - 1, 0] # rank 4 for m in range(cutoff): for n in range(cutoff): for p in range(cutoff): q = p - (m - n) if 0 < q < cutoff: Z[m, n, p, q] = ( R[1, 3] * sqrt[n] / sqrt[q] * Z[m, n - 1, p, q - 1] + R[2, 3] * sqrt[p] / sqrt[q] * Z[m, n, p - 1, q - 1] ) return Z
[docs]@jit(nopython=True) def grad_two_mode_squeezing(T, r, theta): # pragma: no cover """Calculates the gradients of the two-mode squeezing gate with respect to the squeezing magnitude and angle Args: T (array[complex]): array representing the gate r (float): squeezing magnitude theta (float): squeezing angle Returns: tuple[array[complex], array[complex]]: The gradient of the two-mode squeezing gate with respect to r and phi """ cutoff = T.shape[0] dtype = T.dtype sechr = 1.0 / np.cosh(r) tanhr = np.tanh(r) ei = np.exp(1j * theta) eic = np.exp(-1j * theta) sqrt = np.sqrt(np.arange(cutoff, dtype=dtype)) grad_r = np.zeros((cutoff, cutoff, cutoff, cutoff), dtype=dtype) grad_theta = np.zeros((cutoff, cutoff, cutoff, cutoff), dtype=dtype) grad_r[0, 0, 0, 0] = -sechr * tanhr # rank 2 for n in range(1, cutoff): grad_r[n, n, 0, 0] = ( -tanhr * T[n, n, 0, 0] + sqrt[n] * sqrt[n] * ei * sechr**2 * T[n - 1, n - 1, 0, 0] ) grad_theta[n, n, 0, 0] = 1j * ei * tanhr * sqrt[n] * sqrt[n] * T[n - 1, n - 1, 0, 0] # rank 3 for m in range(cutoff): for n in range(m): p = m - n if 0 < p < cutoff: grad_r[m, n, p, 0] = ( -tanhr * T[m, n, p, 0] + sqrt[m] * sqrt[n] * ei * sechr**2 * T[m - 1, n - 1, p, 0] - tanhr * sechr * sqrt[m] * sqrt[p] * T[m - 1, n, p - 1, 0] ) grad_theta[m, n, p, 0] = 1j * ei * tanhr * sqrt[m] * sqrt[n] * T[m - 1, n - 1, p, 0] # rank 4 for m in range(cutoff): for n in range(cutoff): for p in range(cutoff): for q in range(cutoff): grad_r[m, n, p, q] = ( -tanhr * T[m, n, p, q] + sqrt[m] * sqrt[n] * ei * sechr**2 * T[m - 1, n - 1, p, q] - tanhr * sechr * sqrt[m] * sqrt[p] * T[m - 1, n, p - 1, q] - tanhr * sechr * sqrt[n] * sqrt[q] * T[m, n - 1, p, q - 1] - sqrt[p] * sqrt[q] * eic * sechr**2 * T[m, n, p - 1, q - 1] ) grad_theta[m, n, p, q] = ( 1j * ei * tanhr * sqrt[m] * sqrt[n] * T[m - 1, n - 1, p, q] + 1j * eic * tanhr * sqrt[p] * sqrt[q] * T[m, n, p - 1, q - 1] ) return grad_r, grad_theta
[docs]@jit(nopython=True) def beamsplitter(theta, phi, cutoff, dtype=np.complex128): # pragma: no cover r"""Calculates the Fock representation of the beamsplitter. Args: theta (float): transmissivity angle of the beamsplitter. The transmissivity is :math:`t=\cos(\theta)` phi (float): reflection phase of the beamsplitter cutoff (int): Fock ladder cutoff dtype (data type): Specifies the data type used for the calculation Returns: array[float]: The Fock representation of the gate """ sqrt = np.sqrt(np.arange(cutoff, dtype=dtype)) ct = np.cos(theta) st = np.sin(theta) * np.exp(1j * phi) R = np.array( [ [0, 0, ct, -np.conj(st)], [0, 0, st, ct], [ct, st, 0, 0], [-np.conj(st), ct, 0, 0], ] ) Z = np.zeros((cutoff, cutoff, cutoff, cutoff), dtype=dtype) Z[0, 0, 0, 0] = 1.0 # rank 3 for m in range(cutoff): for n in range(cutoff - m): p = m + n if 0 < p < cutoff: Z[m, n, p, 0] = ( R[0, 2] * sqrt[m] / sqrt[p] * Z[m - 1, n, p - 1, 0] + R[1, 2] * sqrt[n] / sqrt[p] * Z[m, n - 1, p - 1, 0] ) # rank 4 for m in range(cutoff): for n in range(cutoff): for p in range(cutoff): q = m + n - p if 0 < q < cutoff: Z[m, n, p, q] = ( R[0, 3] * sqrt[m] / sqrt[q] * Z[m - 1, n, p, q - 1] + R[1, 3] * sqrt[n] / sqrt[q] * Z[m, n - 1, p, q - 1] ) return Z
[docs]@jit(nopython=True) def grad_beamsplitter(T, theta, phi): # pragma: no cover r"""Calculates the gradients of the beamsplitter gate with respect to the transmissivity angle and reflection phase Args: T (array[complex]): array representing the gate theta (float): transmissivity angle of the beamsplitter. The transmissivity is :math:`t=\cos(\theta)` phi (float): reflection phase of the beamsplitter Returns: tuple[array[complex], array[complex]]: The gradient of the beamsplitter gate with respect to theta and phi """ cutoff = T.shape[0] dtype = T.dtype sqrt = np.sqrt(np.arange(cutoff, dtype=dtype)) grad_theta = np.zeros((cutoff, cutoff, cutoff, cutoff), dtype=dtype) grad_phi = np.zeros((cutoff, cutoff, cutoff, cutoff), dtype=dtype) ct = np.cos(theta) st = np.sin(theta) ei = np.exp(1j * phi) eic = np.exp(-1j * phi) # rank 3 for m in range(cutoff): for n in range(cutoff - m): p = m + n if 0 < p < cutoff: grad_theta[m, n, p, 0] = ( -sqrt[m] * sqrt[p] * st * T[m - 1, n, p - 1, 0] + sqrt[n] * sqrt[p] * ei * ct * T[m, n - 1, p - 1, 0] ) grad_phi[m, n, p, 0] = 1j * sqrt[n] * sqrt[p] * ei * st * T[m, n - 1, p - 1, 0] for m in range(cutoff): for n in range(cutoff): for p in range(cutoff): q = m + n - p if 0 < q < cutoff: grad_theta[m, n, p, q] = ( -sqrt[m] * sqrt[p] * st * T[m - 1, n, p - 1, q] - sqrt[n] * sqrt[q] * st * T[m, n - 1, p, q - 1] + sqrt[n] * sqrt[p] * ei * ct * T[m, n - 1, p - 1, q] - sqrt[m] * sqrt[q] * eic * ct * T[m - 1, n, p, q - 1] ) grad_phi[m, n, p, q] = ( 1j * sqrt[n] * sqrt[p] * ei * st * T[m, n - 1, p - 1, q] + 1j * sqrt[m] * sqrt[q] * eic * st * T[m - 1, n, p, q - 1] ) return grad_theta, grad_phi
[docs]@jit(nopython=True) def mzgate(theta, phi, cutoff, dtype=np.complex128): # pragma: no cover r"""Calculates the Fock representation of the Mach-Zehnder interferometer. Args: theta (float): internal phase of the Mach-Zehnder interferometer phi (float): external phase of the Mach-Zehnder interferometer cutoff (int): Fock ladder cutoff dtype (data type): Specifies the data type used for the calculation Returns: array[float]: The Fock representation of the gate """ sqrt = np.sqrt(np.arange(cutoff, dtype=dtype)) v = np.exp(1j * theta) u = np.exp(1j * phi) R = 0.5 * np.array( [ [0, 0, u * (v - 1), 1j * (1 + v)], [0, 0, 1j * u * (1 + v), 1 - v], [u * (v - 1), 1j * u * (1 + v), 0, 0], [1j * (1 + v), 1 - v, 0, 0], ] ) Z = np.zeros((cutoff, cutoff, cutoff, cutoff), dtype=dtype) Z[0, 0, 0, 0] = 1.0 # rank 3 for m in range(cutoff): for n in range(cutoff - m): p = m + n if 0 < p < cutoff: Z[m, n, p, 0] = ( R[0, 2] * sqrt[m] / sqrt[p] * Z[m - 1, n, p - 1, 0] + R[1, 2] * sqrt[n] / sqrt[p] * Z[m, n - 1, p - 1, 0] ) # rank 4 for m in range(cutoff): for n in range(cutoff): for p in range(cutoff): q = m + n - p if 0 < q < cutoff: Z[m, n, p, q] = ( R[0, 3] * sqrt[m] / sqrt[q] * Z[m - 1, n, p, q - 1] + R[1, 3] * sqrt[n] / sqrt[q] * Z[m, n - 1, p, q - 1] ) return Z
[docs]@jit(nopython=True) def grad_mzgate(T, theta, phi): # pragma: no cover r"""Calculates the gradients of the Mach-Zehnder interferometer with respect to the transmissivity angle and reflection phase Args: T (array[complex]): array representing the gate theta (float): internal of the mzgate phi (float): external phase of the mzgate Returns: tuple[array[complex], array[complex]]: The gradient of the mzgate gate with respect to theta and phi """ cutoff = T.shape[0] dtype = T.dtype sqrt = np.sqrt(np.arange(cutoff, dtype=dtype)) grad_theta = np.zeros((cutoff, cutoff, cutoff, cutoff), dtype=dtype) grad_phi = np.zeros((cutoff, cutoff, cutoff, cutoff), dtype=dtype) v = np.exp(1j * theta) u = np.exp(1j * phi) # rank 3 for m in range(cutoff): for n in range(cutoff - m): p = m + n if 0 < p < cutoff: grad_theta[m, n, p, 0] = ( 1j * 0.5 * u * v * sqrt[m] * sqrt[p] * T[m - 1, n, p - 1, 0] - 0.5 * u * v * sqrt[n] * sqrt[p] * T[m, n - 1, p - 1, 0] ) grad_phi[m, n, p, 0] = (1j * 0.5 * u * v - 1j * 0.5 * u) * sqrt[m] * sqrt[p] * T[ m - 1, n, p - 1, 0 ] - (0.5 * u + 0.5 * u * v) * sqrt[n] * sqrt[p] * T[m, n - 1, p - 1, 0] for m in range(cutoff): for n in range(cutoff): for p in range(cutoff): q = m + n - p if 0 < q < cutoff: grad_theta[m, n, p, q] = ( 1j * 0.5 * u * v * sqrt[m] * sqrt[p] * T[m - 1, n, p - 1, q] - 0.5 * u * v * sqrt[n] * sqrt[p] * T[m, n - 1, p - 1, q] - 0.5 * v * sqrt[m] * sqrt[q] * T[m - 1, n, p, q - 1] - 1j * 0.5 * v * sqrt[n] * sqrt[q] * T[m, n - 1, p, q - 1] ) grad_phi[m, n, p, q] = (1j * 0.5 * u * v - 1j * 0.5 * u) * sqrt[m] * sqrt[ p ] * T[m - 1, n, p - 1, q] - (0.5 * u + 0.5 * u * v) * sqrt[n] * sqrt[p] * T[ m, n - 1, p - 1, q ] return grad_theta, grad_phi