Source code for mrmustard.physics.gaussian
# Copyright 2021 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.
"""
This module contains functions for performing calculations on objects in the Gaussian representations.
"""
from typing import Any, Optional, Sequence, Tuple, Union
from mrmustard import math, settings
from mrmustard.math.tensor_wrappers.xptensor import XPMatrix, XPVector
from mrmustard.utils.typing import Matrix, Scalar, Vector
# ~~~~~~
# States
# ~~~~~~
[docs]
def vacuum_cov(num_modes: int) -> Matrix:
r"""Returns the real covariance matrix of the vacuum state.
Args:
num_modes (int): number of modes
Returns:
Matrix: vacuum covariance matrix
"""
return math.eye(num_modes * 2, dtype=math.float64) * settings.HBAR / 2
[docs]
def vacuum_means(num_modes: int) -> Tuple[Matrix, Vector]:
r"""Returns the real covariance matrix and real means vector of the vacuum state.
Args:
num_modes (int): number of modes
Returns:
Matrix, Vector: thermal state covariance matrix or means vector
"""
return displacement(
math.zeros(num_modes, dtype="float64"),
math.zeros(num_modes, dtype="float64"),
)
[docs]
def squeezed_vacuum_cov(r: Vector, phi: Vector) -> Matrix:
r"""Returns the real covariance matrix and real means vector of a squeezed vacuum state.
The dimension depends on the dimensions of ``r`` and ``phi``.
Args:
r (vector): squeezing magnitude
phi (vector): squeezing angle
Returns:
Matrix, Vector: thermal state covariance matrix or means vector
"""
S = squeezing_symplectic(r, phi)
return math.matmul(S, math.transpose(S)) * settings.HBAR / 2
[docs]
def thermal_cov(nbar: Vector) -> Tuple[Matrix, Vector]:
r"""Returns the real covariance matrix and real means vector of a thermal state.
The dimension depends on the dimensions of ``nbar``.
Args:
nbar (vector): average number of photons per mode
Returns:
Matrix, Vector: thermal state covariance matrix or means vector
"""
g = (2 * math.atleast_1d(nbar) + 1) * settings.HBAR / 2
return math.diag(math.concat([g, g], axis=-1))
[docs]
def two_mode_squeezed_vacuum_cov(r: Vector, phi: Vector) -> Matrix:
r"""Returns the real covariance matrix and real means vector of a two-mode squeezed vacuum state.
The dimension depends on the dimensions of ``r`` and ``phi``.
Args:
r (vector): squeezing magnitude
phi (vector): squeezing angle
Returns:
Matrix: two-mode squeezed state covariance matrix
Vector: two-mode squeezed state means vector
"""
S = two_mode_squeezing_symplectic(r, phi)
return math.matmul(S, math.transpose(S)) * settings.HBAR / 2
[docs]
def gaussian_cov(symplectic: Matrix, eigenvalues: Vector = None) -> Matrix:
r"""Returns the covariance matrix of a Gaussian state.
Args:
symplectic (Tensor): symplectic matrix of a channel
eigenvalues (vector): symplectic eigenvalues
Returns:
Tensor: covariance matrix of the Gaussian state
"""
if eigenvalues is None:
return math.matmul(symplectic, math.transpose(symplectic))
return math.matmul(
math.matmul(symplectic, math.diag(math.concat([eigenvalues, eigenvalues], axis=0))),
math.transpose(symplectic),
)
# ~~~~~~~~~~~~~~~~~~~~~~~~
# Unitary transformations
# ~~~~~~~~~~~~~~~~~~~~~~~~
[docs]
def rotation_symplectic(angle: Union[Scalar, Vector]) -> Matrix:
r"""Symplectic matrix of a rotation gate.
The dimension depends on the dimension of the angle.
Args:
angle (scalar or vector): rotation angles
Returns:
Tensor: symplectic matrix of a rotation gate
"""
angle = math.atleast_1d(angle)
num_modes = angle.shape[-1]
x = math.cos(angle)
y = math.sin(angle)
return (
math.diag(math.concat([x, x], axis=0))
+ math.diag(-y, k=num_modes)
+ math.diag(y, k=-num_modes)
)
[docs]
def squeezing_symplectic(r: Union[Scalar, Vector], phi: Union[Scalar, Vector]) -> Matrix:
r"""Symplectic matrix of a squeezing gate.
The dimension depends on the dimension of ``r`` and ``phi``.
Args:
r (scalar or vector): squeezing magnitude
phi (scalar or vector): rotation parameter
Returns:
Tensor: symplectic matrix of a squeezing gate
"""
r = math.atleast_1d(r, math.float64)
phi = math.atleast_1d(phi, math.float64)
if r.shape[-1] == 1:
r = math.tile(r, phi.shape)
if phi.shape[-1] == 1:
phi = math.tile(phi, r.shape)
num_modes = phi.shape[-1]
cp = math.cos(phi)
sp = math.sin(phi)
ch = math.cosh(r)
sh = math.sinh(r)
cpsh = cp * sh
spsh = sp * sh
return (
math.diag(math.concat([ch - cpsh, ch + cpsh], axis=0))
+ math.diag(-spsh, k=num_modes)
+ math.diag(-spsh, k=-num_modes)
)
[docs]
def displacement(x: Union[Scalar, Vector], y: Union[Scalar, Vector]) -> Vector:
r"""Returns the displacement vector for a displacement by :math:`alpha = x + iy`.
The dimension depends on the dimensions of ``x`` and ``y``.
Args:
x (scalar or vector): real part of displacement (in units of :math:`\sqrt{\hbar}`)
y (scalar or vector): imaginary part of displacement (in units of :math:`\sqrt{\hbar}`)
Returns:
Vector: displacement vector of a displacement gate
"""
x = math.atleast_1d(x, math.float64)
y = math.atleast_1d(y, math.float64)
if x.shape[-1] == 1:
x = math.tile(x, y.shape)
if y.shape[-1] == 1:
y = math.tile(y, x.shape)
return math.sqrt(2 * settings.HBAR, dtype=x.dtype) * math.concat([x, y], axis=0)
[docs]
def beam_splitter_symplectic(theta: Scalar, phi: Scalar) -> Matrix:
r"""Symplectic matrix of a Beam-splitter gate.
The dimension is :math:`4\times 4`.
Args:
theta: transmissivity parameter
phi: phase parameter
Returns:
Matrix: symplectic (orthogonal) matrix of a beam-splitter gate
"""
ct = math.cos(theta)
st = math.sin(theta)
cp = math.cos(phi)
sp = math.sin(phi)
zero = math.zeros_like(theta)
return math.astensor(
[
[ct, -cp * st, zero, -sp * st],
[cp * st, ct, -sp * st, zero],
[zero, sp * st, ct, -cp * st],
[sp * st, zero, cp * st, ct],
]
)
[docs]
def mz_symplectic(phi_a: Scalar, phi_b: Scalar, internal: bool = False) -> Matrix:
r"""Symplectic matrix of a Mach-Zehnder gate.
It supports two conventions:
* if ``internal=True``, both phases act inside the interferometer:
``phi_a`` on the upper arm, ``phi_b`` on the lower arm;
* if `internal = False` (default), both phases act on the upper arm:
``phi_a`` before the first BS, ``phi_b`` after the first BS.
Args:
phi_a (float): first phase
phi_b (float): second phase
internal (bool): whether phases are in the internal arms (default is False)
Returns:
Matrix: symplectic (orthogonal) matrix of a Mach-Zehnder interferometer
"""
ca = math.cos(phi_a)
sa = math.sin(phi_a)
cb = math.cos(phi_b)
sb = math.sin(phi_b)
cp = math.cos(phi_a + phi_b)
sp = math.sin(phi_a + phi_b)
if internal:
return 0.5 * math.astensor(
[
[ca - cb, -sa - sb, sb - sa, -ca - cb],
[-sa - sb, cb - ca, -ca - cb, sa - sb],
[sa - sb, ca + cb, ca - cb, -sa - sb],
[ca + cb, sb - sa, -sa - sb, cb - ca],
]
)
return 0.5 * math.astensor(
[
[cp - ca, -sb, sa - sp, -1 - cb],
[-sa - sp, 1 - cb, -ca - cp, sb],
[sp - sa, 1 + cb, cp - ca, -sb],
[cp + ca, -sb, -sa - sp, 1 - cb],
]
)
[docs]
def two_mode_squeezing_symplectic(r: Scalar, phi: Scalar) -> Matrix:
r"""Symplectic matrix of a two-mode squeezing gate.
The dimension is :math:`4\times 4`.
Args:
r (float): squeezing magnitude
phi (float): rotation parameter
Returns:
Matrix: symplectic matrix of a two-mode squeezing gate
"""
cp = math.cast(math.cos(phi), math.float64)
sp = math.cast(math.sin(phi), math.float64)
ch = math.cast(math.cosh(r), math.float64)
sh = math.cast(math.sinh(r), math.float64)
zero = math.cast(math.zeros_like(math.asnumpy(r)), math.float64)
return math.astensor(
[
[ch, cp * sh, zero, sp * sh],
[cp * sh, ch, sp * sh, zero],
[zero, sp * sh, ch, -cp * sh],
[sp * sh, zero, -cp * sh, ch],
]
)
[docs]
def quadratic_phase(s: Scalar):
r"""Quadratic phase single mode gate.
.. math::
P = \exp(i s q^2 / 2 \hbar)
Reference: https://strawberryfields.ai/photonics/conventions/gates.html
Args:
s (float): interaction strength
Returns:
Tensor: the :math:`P(s)` matrix (in ``xxpp`` ordering)
"""
return math.astensor(
[
[1, 0],
[s, 1],
]
)
[docs]
def controlled_Z(g: Scalar):
r"""Controlled PHASE gate of two-gaussian modes.
.. math::
C_Z = \exp(ig q_1 \otimes q_2 / \hbar).
Reference: https://arxiv.org/pdf/2110.03247.pdf, Equation 8.
https://arxiv.org/pdf/1110.3234.pdf, Equation 161.
Args:
g (float): interaction strength
Returns:
the C_Z(g) matrix (in xxpp ordering)
"""
return math.astensor(
[
[1, 0, 0, 0],
[0, 1, 0, 0],
[0, g, 1, 0],
[g, 0, 0, 1],
]
)
[docs]
def controlled_X(g: Scalar):
r"""Controlled NOT gate of two-gaussian modes.
.. math::
C_X = \exp(ig q_1 \otimes p_2).
Reference: https://arxiv.org/pdf/2110.03247.pdf, Equation 9.
Args:
g (float): interaction strength
Returns:
the C_X(g) matrix (in xxpp ordering)
"""
return math.astensor(
[
[1, 0, 0, 0],
[g, 1, 0, 0],
[0, 0, 1, -g],
[0, 0, 0, 1],
]
)
# ~~~~~~~~~~~~~
# CPTP channels
# ~~~~~~~~~~~~~
[docs]
def CPTP(
cov: Matrix,
means: Vector,
X: Matrix,
Y: Matrix,
d: Vector,
state_modes: Sequence[int],
transf_modes: Sequence[int],
) -> Tuple[Matrix, Vector]:
r"""Returns the cov matrix and means vector of a state after undergoing a CPTP channel.
Computed as ``cov = X \cdot cov \cdot X^T + Y`` and ``d = X \cdot means + d``.
If the channel is single-mode, ``modes`` can contain ``M`` modes to apply the channel to,
otherwise it must contain as many modes as the number of modes in the channel.
Args:
cov (Matrix): covariance matrix
means (Vector): means vector
X (Matrix): the X matrix of the CPTP channel
Y (Matrix): noise matrix of the CPTP channel
d (Vector): displacement vector of the CPTP channel
state_modes (Sequence[int]): modes the state is defined on
transf_modes (Sequence[int]): modes on which the channel acts
Returns:
Tuple[Matrix, Vector]: the covariance matrix and the means vector of the state after the CPTP channel
"""
if not set(transf_modes).issubset(state_modes):
raise ValueError(
f"The channel should act on a subset of the state modes ({transf_modes} is not a subset of {state_modes})"
)
# if single-mode channel, apply to all modes indicated in `modes`
if X is not None and X.shape[-1] == 2:
X = math.single_mode_to_multimode_mat(X, len(transf_modes))
if Y is not None and Y.shape[-1] == 2:
Y = math.single_mode_to_multimode_mat(Y, len(transf_modes))
if d is not None and d.shape[-1] == 2:
d = math.single_mode_to_multimode_vec(d, len(transf_modes))
indices = [
state_modes.index(i) for i in transf_modes
] # TODO: do this when calling the method instead of here?
cov = math.left_matmul_at_modes(X, cov, indices)
cov = math.right_matmul_at_modes(cov, math.transpose(X), indices)
cov = math.add_at_modes(cov, Y, indices)
means = math.matvec_at_modes(X, means, indices)
means = math.add_at_modes(means, d, indices)
return cov, means
[docs]
def loss_XYd(
transmissivity: Union[Scalar, Vector], nbar: Union[Scalar, Vector]
) -> Tuple[Matrix, Matrix, None]:
r"""Returns the ``X``, ``Y`` matrices and the ``d`` vector for the noisy loss (attenuator) channel.
.. math::
X = math.sqrt(transmissivity)
Y = (1-transmissivity) * (2 * nbar + 1) * hbar / 2
Reference: Alessio Serafini - Quantum Continuous Variables (5.77, p. 108)
Args:
transmissivity (float): value of the transmissivity, must be between 0 and 1
nbar (float): photon number expectation value in the environment (0 for pure loss channel)
Returns:
Tuple[Matrix, Matrix, None]: the ``X``, ``Y`` matrices and the ``d`` vector for the noisy
loss channel
"""
if math.any(transmissivity < 0) or math.any(transmissivity > 1):
raise ValueError("transmissivity must be between 0 and 1")
x = math.sqrt(transmissivity)
X = math.diag(math.concat([x, x], axis=0))
y = (1 - transmissivity) * (2 * nbar + 1) * settings.HBAR / 2
Y = math.diag(math.concat([y, y], axis=0))
return X, Y, None
[docs]
def amp_XYd(gain: Union[Scalar, Vector], nbar: Union[Scalar, Vector]) -> Matrix:
r"""Returns the ``X``, ``Y`` matrices and the d vector for the noisy amplifier channel.
.. math::
X = math.sqrt(gain)
Y = (gain-1) * (2 * nbar + 1) * hbar / 2
Reference: Alessio Serafini - Quantum Continuous Variables (5.77, p. 111)
The quantum limited amplifier channel is recovered for ``nbar = 0.0``.
Args:
gain (float): value of the gain > 1
nbar (float): photon number expectation value in the environment (0 for quantum
limited amplifier)
Returns:
Tuple[Matrix, Vector]: the ``X``, ``Y`` matrices and the ``d`` vector for the noisy
amplifier channel.
"""
if math.any(gain < 1):
raise ValueError("Gain must be larger than 1")
x = math.sqrt(gain)
X = math.diag(math.concat([x, x], axis=0))
y = (gain - 1) * (2 * nbar + 1) * settings.HBAR / 2
Y = math.diag(math.concat([y, y], axis=0))
return X, Y, None
[docs]
def noise_Y(noise: Union[Scalar, Vector]) -> Matrix:
r"""Returns the ``X``, ``Y`` matrices and the d vector for the additive noise channel ``(Y = noise * (\hbar / 2) * I)``
Args:
noise (float): number of photons in the thermal state
Returns:
Tuple[None, Matrix, None]: the ``X``, ``Y`` matrices and the ``d`` vector of the noise channel.
"""
return math.diag(math.concat([noise, noise], axis=0)) * settings.HBAR / 2
[docs]
def compose_channels_XYd(
X1: Matrix, Y1: Matrix, d1: Vector, X2: Matrix, Y2: Matrix, d2: Vector
) -> Tuple[Matrix, Matrix, Vector]:
r"""Returns the combined ``X``, ``Y``, and ``d`` for two CPTP channels.
Args:
X1 (Matrix): the ``X`` matrix of the first CPTP channel
Y1 (Matrix): the ``Y`` matrix of the first CPTP channel
d1 (Vector): the displacement vector of the first CPTP channel
X2 (Matrix): the ``X`` matrix of the second CPTP channel
Y2 (Matrix): the ``Y`` matrix of the second CPTP channel
d2 (Vector): the displacement vector of the second CPTP channel
Returns:
Tuple[Matrix, Matrix, Vector]: the combined ``X``, ``Y``, and ``d`` matrices
"""
if X1 is None:
X = X2
elif X2 is None:
X = X1
else:
X = math.matmul(X2, X1)
if Y1 is None:
Y = Y2
elif Y2 is None:
Y = Y1
else:
Y = math.matmul(math.matmul(X2, Y1), X2) + Y2
if d1 is None:
d = d2
elif d2 is None:
d = d1
else:
d = math.matmul(X2, d1) + d2
return X, Y, d
# ~~~~~~~~~~~~~~~
# non-TP channels
# ~~~~~~~~~~~~~~~
[docs]
def general_dyne(
cov: Matrix,
means: Vector,
proj_cov: Matrix,
proj_means: Optional[Vector] = None,
modes: Optional[Sequence[int]] = None,
) -> Tuple[Scalar, Matrix, Vector]:
r"""Returns the results of a general-dyne measurement. If ``proj_means`` are not provided
(as ``None``), they are sampled from the probability distribution.
Args:
cov (Matrix): covariance matrix of the state being measured [units of `2\hbar`]
means (Vector): means vector of the state being measured [units of `\sqrt(\hbar)`]
proj_cov (Matrix): covariance matrix of the state being projected onto [units of `2\hbar`]
proj_means (Optional Vector): means vector of the state being projected onto
(i.e. the measurement outcome) [units of `\sqrt(\hbar)`]. If not provided, the means vector
is sampled from the generaldyne probability distribution.
modes (Optional Sequence[int]): modes being measured (modes are indexed from 0 to num_modes-1),
if modes are not provided then the first modes (according to the size of ``cov``) are measured.
Returns:
Tuple[Scalar, Scalar, Matrix, Vector]:
outcome (sampled means vector of the measured subsystem) [units of `\sqrt(\hbar)`],
oucome probability [units of `\hbar**N`],
post-measurement covariace [units of `2\hbar`]
post-measurement means vector [units of `\sqrt{\hbar}`].
"""
N, M = cov.shape[-1] // 2, proj_cov.shape[-1] // 2
# Bmodes are the modes being measured and Amodes are the leftover modes
Bmodes = modes or list(range(M))
Amodes = list(set(list(range(N))) - set(Bmodes))
A, B, AB = partition_cov(cov, Amodes)
a, b = partition_means(means, Amodes)
reduced_cov = B + proj_cov
# covariances are divided by 2 to match tensorflow and MrMustard conventions
# (MrMustard uses Serafini convention where `sigma_MM = 2 sigma_TF`)
if proj_means is None:
pdf = math.MultivariateNormalTriL(loc=b, scale_tril=math.cholesky(reduced_cov / 2))
outcome = (
pdf.sample(dtype=cov.dtype) if proj_means is None else math.cast(proj_means, cov.dtype)
)
prob = pdf.prob(outcome)
else:
# If the projector is already given: proj_means
# use the formula 5.139 in Serafini - Quantum Continuous Variables
# fixed by -0.5 on the exponential, added hbar and removed pi due to different convention
outcome = proj_means
prob = (
settings.HBAR**M
* math.exp(
-0.5 * math.sum(math.solve(reduced_cov, (proj_means - b)) * (proj_means - b))
)
/ math.sqrt(math.det(reduced_cov))
)
# calculate conditional output state of unmeasured modes
num_remaining_modes = N - M
if num_remaining_modes == 0:
return outcome, prob, None, None
AB_inv = math.matmul(AB, math.inv(reduced_cov))
new_cov = A - math.matmul(AB_inv, math.transpose(AB))
new_means = a + math.matvec(AB_inv, outcome - b)
return outcome, prob, new_cov, new_means
# ~~~~~~~~~
# utilities
# ~~~~~~~~~
[docs]
def number_means(cov: Matrix, means: Vector) -> Vector:
r"""Returns the photon number means vector given a Wigner covariance matrix and a means vector.
Args:
cov: the Wigner covariance matrix
means: the Wigner means vector
Returns:
Vector: the photon number means vector
"""
N = means.shape[-1] // 2
return (
means[:N] ** 2
+ means[N:] ** 2
+ math.diag_part(cov[:N, :N])
+ math.diag_part(cov[N:, N:])
- settings.HBAR
) / (2 * settings.HBAR)
[docs]
def number_cov(cov: Matrix, means: Vector) -> Matrix:
r"""Returns the photon number covariance matrix given a Wigner covariance matrix and a means vector.
Args:
cov: the Wigner covariance matrix
means: the Wigner means vector
Returns:
Matrix: the photon number covariance matrix
"""
N = means.shape[-1] // 2
mCm = cov * means[:, None] * means[None, :]
dd = math.diag(math.diag_part(mCm[:N, :N] + mCm[N:, N:] + mCm[:N, N:] + mCm[N:, :N])) / (
2 * settings.HBAR**2 # TODO: sum(diag_part) is better than diag_part(sum)
)
CC = (cov**2 + mCm) / (2 * settings.HBAR**2)
return (
CC[:N, :N] + CC[N:, N:] + CC[:N, N:] + CC[N:, :N] + dd - 0.25 * math.eye(N, dtype=CC.dtype)
)
[docs]
def trace(cov: Matrix, means: Vector, Bmodes: Sequence[int]) -> Tuple[Matrix, Vector]:
r"""Returns the covariances and means after discarding the specified modes.
Args:
cov (Matrix): covariance matrix
means (Vector): means vector
Bmodes (Sequence[int]): modes to discard
Returns:
Tuple[Matrix, Vector]: the covariance matrix and the means vector after discarding the specified modes
"""
N = len(cov) // 2
Aindices = math.astensor(
[i for i in range(N) if i not in Bmodes] + [i + N for i in range(N) if i not in Bmodes],
dtype=math.int32,
)
A_cov_block = math.gather(math.gather(cov, Aindices, axis=0), Aindices, axis=1)
A_means_vec = math.gather(means, Aindices)
return A_cov_block, A_means_vec
[docs]
def partition_cov(cov: Matrix, Amodes: Sequence[int]) -> Tuple[Matrix, Matrix, Matrix]:
r"""Partitions the covariance matrix into the ``A`` and ``B`` subsystems and the AB coherence block.
Args:
cov (Matrix): the covariance matrix
Amodes (Sequence[int]): the modes of system ``A``
Returns:
Tuple[Matrix, Matrix, Matrix]: the cov of ``A``, the cov of ``B`` and the AB block
"""
N = cov.shape[-1] // 2
Bindices = math.cast(
[i for i in range(N) if i not in Amodes] + [i + N for i in range(N) if i not in Amodes],
"int32",
)
Aindices = math.cast(Amodes + [i + N for i in Amodes], "int32")
A_block = math.gather(math.gather(cov, Aindices, axis=1), Aindices, axis=0)
B_block = math.gather(math.gather(cov, Bindices, axis=1), Bindices, axis=0)
AB_block = math.gather(math.gather(cov, Bindices, axis=1), Aindices, axis=0)
return A_block, B_block, AB_block
[docs]
def partition_means(means: Vector, Amodes: Sequence[int]) -> Tuple[Vector, Vector]:
r"""Partitions the means vector into the ``A`` and ``B`` subsystems.
Args:
means (Vector): the means vector
Amodes (Sequence[int]): the modes of system ``A``
Returns:
Tuple[Vector, Vector]: the means of ``A`` and the means of ``B``
"""
N = len(means) // 2
Bindices = math.cast(
[i for i in range(N) if i not in Amodes] + [i + N for i in range(N) if i not in Amodes],
"int32",
)
Aindices = math.cast(Amodes + [i + N for i in Amodes], "int32")
return math.gather(means, Aindices), math.gather(means, Bindices)
[docs]
def purity(cov: Matrix) -> Scalar:
r"""Returns the purity of the state with the given covariance matrix.
Args:
cov (Matrix): the covariance matrix
Returns:
float: the purity
"""
return 1 / math.sqrt(math.det((2 / settings.HBAR) * cov))
[docs]
def symplectic_eigenvals(cov: Matrix) -> Any:
r"""Returns the sympletic eigenspectrum of a covariance matrix.
For a pure state, we expect the sympletic eigenvalues to be 1.
Args:
cov (Matrix): the covariance matrix
Returns:
List[float]: the sympletic eigenvalues
"""
J = math.J(cov.shape[-1] // 2) # create a sympletic form
M = math.matmul(1j * J, cov * (2 / settings.HBAR))
vals = math.abs(math.eigvals(M)) # compute the eigenspectrum
vals = math.sort(vals)
return vals[::2] # return the even eigenvalues
[docs]
def von_neumann_entropy(cov: Matrix) -> float:
r"""Returns the Von Neumann entropy.
For a pure state, we expect the Von Neumann entropy to be 0.
Reference: (https://arxiv.org/pdf/1110.3234.pdf), Equations 46-47.
Args:
cov (Matrix): the covariance matrix
Returns:
float: the Von Neumann entropy
"""
def g(x):
return math.xlogy((x + 1) / 2, (x + 1) / 2) - math.xlogy((x - 1) / 2, (x - 1) / 2 + 1e-9)
symp_vals = symplectic_eigenvals(cov)
entropy = math.sum(g(symp_vals))
return entropy
[docs]
def fidelity(mu1: Vector, cov1: Matrix, mu2: Vector, cov2: Matrix) -> float:
r"""Returns the fidelity of two gaussian states.
Reference: `arXiv:2102.05748 <https://arxiv.org/pdf/2102.05748.pdf>`_, equations 95-99.
Note that we compute the square of equation 98.
Args:
mu1 (Vector): the means vector of state 1
mu2 (Vector): the means vector of state 2
cov1 (Matrix): the covariance matrix of state 1
cov1 (Matrix): the covariance matrix of state 2
Returns:
float: the fidelity
"""
cov1 = math.cast(cov1 / settings.HBAR, "complex128") # convert to units where hbar = 1
cov2 = math.cast(cov2 / settings.HBAR, "complex128") # convert to units where hbar = 1
mu1 = math.cast(mu1, "complex128")
mu2 = math.cast(mu2, "complex128")
deltar = (mu2 - mu1) / math.sqrt(
settings.HBAR, dtype=mu1.dtype
) # convert to units where hbar = 1
J = math.J(cov1.shape[0] // 2)
I = math.eye(cov1.shape[0])
J = math.cast(J, "complex128")
I = math.cast(I, "complex128")
cov12_inv = math.inv(cov1 + cov2)
V = math.transpose(J) @ cov12_inv @ ((1 / 4) * J + cov2 @ J @ cov1)
W = -2 * (V @ (1j * J))
W_inv = math.inv(W)
matsqrtm = math.sqrtm(
I - W_inv @ W_inv
) # this also handles the case where the input matrix is close to zero
f0_top = math.det((matsqrtm + I) @ (W @ (1j * J)))
f0_bot = math.det(cov1 + cov2)
f0 = math.sqrt(f0_top / f0_bot) # square of equation 98
dot = math.sum(
math.transpose(deltar) * math.matvec(cov12_inv, deltar)
) # computing (mu2-mu1)/sqrt(hbar).T @ cov12_inv @ (mu2-mu1)/sqrt(hbar)
_fidelity = f0 * math.exp((-1 / 2) * dot) # square of equation 95
return math.real(_fidelity)
[docs]
def physical_partial_transpose(cov: Matrix, modes: Sequence[int]) -> Matrix:
r"""Returns the covariance matrix that corresponds to applying the partial
transposition on the density matrix of a given set of modes.
Reference: `https://arxiv.org/abs/quant-ph/9909044 <https://arxiv.org/abs/quant-ph/9909044>`_, Equation 1, 5.
Args:
cov (Matrix): the covariance matrix
modes (Sequence[int]): the modes of system on which transposition is applied
Returns:
Matrix: the covariance matrix corresponding to the partially transposed state
"""
m, _ = cov.shape
num_modes = m // 2
mat = [1.0] * m
for i in modes:
mat[i + num_modes] = -1.0
mat = math.astensor(mat, dtype="float64")
return cov * mat[:, None] * mat[None, :]
[docs]
def log_negativity(cov: Matrix) -> float:
r"""Returns the log_negativity of a Gaussian state.
Reference: `https://arxiv.org/abs/quant-ph/0102117 <https://arxiv.org/abs/quant-ph/0102117>`_ , Equation 57, 61.
Args:
cov (Matrix): the covariance matrix
Returns:
float: the log-negativity
"""
vals = symplectic_eigenvals(cov) / (settings.HBAR / 2)
vals_filtered = math.boolean_mask(
vals, vals < 1.0
) # Get rid of terms that would lead to zero contribution.
if len(vals_filtered) > 0:
return -math.sum(
math.log(vals_filtered) / math.cast(math.log(2.0), dtype=vals_filtered.dtype)
)
return 0
[docs]
def join_covs(covs: Sequence[Matrix]) -> Tuple[Matrix, Vector]:
r"""Joins the given covariance matrices into a single covariance matrix.
Args:
covs (Sequence[Matrix]): the covariance matrices
Returns:
Matrix: the joined covariance matrix
"""
modes = list(range(len(covs[0]) // 2))
cov = XPMatrix.from_xxpp(covs[0], modes=(modes, modes), like_1=True)
for _, c in enumerate(covs[1:]):
modes = list(range(cov.num_modes, cov.num_modes + c.shape[-1] // 2))
cov = cov @ XPMatrix.from_xxpp(c, modes=(modes, modes), like_1=True)
return cov.to_xxpp()
[docs]
def join_means(means: Sequence[Vector]) -> Vector:
r"""Joins the given means vectors into a single means vector.
Args:
means (Sequence[Vector]): the means vectors
Returns:
Vector: the joined means vector
"""
mean = XPVector.from_xxpp(means[0], modes=list(range(len(means[0]) // 2)))
for _, m in enumerate(means[1:]):
mean = mean + XPVector.from_xxpp(
m, modes=list(range(mean.num_modes, mean.num_modes + len(m) // 2))
)
return mean.to_xxpp()
[docs]
def symplectic_inverse(S: Matrix) -> Matrix:
r"""Returns the inverse of a symplectic matrix.
Args:
S (Matrix): the symplectic matrix
Returns:
Matrix: the inverse of the symplectic matrix
"""
S = math.reshape(S, (S.shape[0] // 2, 2, S.shape[1] // 2, 2))
S = math.transpose(S, (1, 3, 0, 2))
return math.block(
[
[math.transpose(S[1, 1]), -math.transpose(S[0, 1])],
[-math.transpose(S[1, 0]), math.transpose(S[0, 0])],
]
)
[docs]
def XYd_dual(X: Matrix, Y: Matrix, d: Vector):
r"""Returns the dual channel ``(X, Y, d)``.
Args:
X (Matrix): the ``X`` matrix
Y (Matrix): the ``Y`` noise matrix
d (Vector): the displacement vector
Returns:
Tuple[Matrix, Matrix, Vector]: ``(X_dual, Y_dual, d_dual)``
"""
X_dual = math.inv(X) if X is not None else None
Y_dual = Y
d_dual = d
if Y is not None:
Y_dual = (
math.matmul(X_dual, math.matmul(Y, math.transpose(X_dual))) if X_dual is not None else Y
)
if d is not None:
d_dual = math.matvec(X_dual, d) if X_dual is not None else d
return X_dual, Y_dual, d_dual
_modules/mrmustard/physics/gaussian
Download Python script
Download Notebook
View on GitHub