Utilities

Module name: strawberryfields.utils

This module defines and implements several utility functions and language extensions that complement StrawberryFields.

Classical processing functions

These functions provide common mathematical operations that may be required for classical processing of measured modes input to other gates. They may be used as follows:

MeasureX | q[0]
Xgate(scale(q[0], sqrt(0.5))) | q[1]

Available classical processing functions include:

neg(x) Negates a measured value.
mag(x) Returns the magnitude \(|z|\) of a measured value.
phase(x) Returns the phase \(\phi\) of a measured value \(z=re^{i\phi}\).
scale(x, a) Scales the measured value by factor a.
shift(x, b) Shifts the measured value by factor b.
scale_shift(x, a, b) Scales the measured value by factor a then shifts the result by b.
power(x, a) Raises the measured value to power a.

If more advanced classical processing is required, custom classical processing functions can be created using the strawberryfields.convert() decorator.

NumPy state functions

These functions allow the calculation of various quantum states in either the Fock basis (a one-dimensional array indexed by Fock state) or the Gaussian basis (returning a vector of means and covariance matrix). These state calculations are NOT done in the simulators, but rather in NumPy.

These are useful for generating states for use in calculating the fidelity of simulations.

squeezed_cov(r, phi[, hbar]) Returns the squeezed covariance matrix of a squeezed state
vacuum_state([basis, fock_dim, hbar]) Returns the vacuum state
coherent_state(a[, basis, fock_dim, hbar]) Returns the coherent state
squeezed_state(r, p[, basis, fock_dim, hbar]) Returns the squeezed state
displaced_squeezed_state(a, r, phi[, basis, …]) Returns the squeezed coherent state
fock_state(n[, fock_dim]) Returns the Fock state
cat_state(a[, p, fock_dim]) Returns the cat state

Random functions

These functions generate random numbers and matrices corresponding to various quantum states and operations.

randnc(*arg) Normally distributed array of random complex numbers.
random_covariance(N[, hbar, pure, block_diag]) Random covariance matrix.
random_symplectic(N[, passive, block_diag]) Random symplectic matrix representing a Gaussian transformation.
random_interferometer(N[, real]) Random unitary matrix representing an interferometer.

Decorators

The operation decorator allows functions containing quantum operations acting on a qumode to be used as an operation itself within a Program context.

operation(ns) Groups a sequence of gates into a single operation to be used within a Program context.

Program functions

These functions act on Program instances, returning or extracting information from the quantum circuit.

For example, these might be used as follows:

prog = sf.Program(2)
with prog.context as q:
    BSgate(0.543, 0.123) | (q[0], q[1])

U = extract_unitary(prog, cutoff_dim=10)

In this example, U is a unitary array representing the quantum circuit prog in the Fock basis (here, a single beamsplitter).

is_unitary(prog) True iff all the operations in the program are unitary.
is_channel(prog) True iff all the operations in the program can be represented as quantum channels.
extract_unitary(prog, cutoff_dim, …) Numerical array representation of a unitary quantum circuit.
extract_channel(prog, cutoff_dim, …) Numerical array representation of the channel corresponding to a quantum circuit.

Code details

strawberryfields.utils.neg(x)[source]

Negates a measured value.

Parameters:x (RegRef) – mode that has been previously measured
strawberryfields.utils.mag(x)[source]

Returns the magnitude \(|z|\) of a measured value.

Parameters:x (RegRef) – mode that has been previously measured
strawberryfields.utils.phase(x)[source]

Returns the phase \(\phi\) of a measured value \(z=re^{i\phi}\).

Parameters:x (RegRef) – mode that has been previously measured
strawberryfields.utils.scale(x, a)[source]

Scales the measured value by factor a.

Parameters:
  • x (RegRef) – mode that has been previously measured
  • a (float) – scaling factor
strawberryfields.utils.shift(x, b)[source]

Shifts the measured value by factor b.

Parameters:
  • x (RegRef) – mode that has been previously measured
  • b (float) – shifting factor
strawberryfields.utils.scale_shift(x, a, b)[source]

Scales the measured value by factor a then shifts the result by b.

\[u' = au + b\]
Parameters:
  • x (RegRef) – mode that has been previously measured
  • a (float) – scaling factor
  • b (float) – shifting factor
strawberryfields.utils.power(x, a)[source]

Raises the measured value to power a.

Parameters:
  • x (RegRef) – mode that has been previously measured
  • a (float) – the exponent of x; note that a can be negative and fractional
strawberryfields.utils.squeezed_cov(r, phi, hbar=2)[source]

Returns the squeezed covariance matrix of a squeezed state

Parameters:
  • r (complex) – the squeezing magnitude
  • p (float) – the squeezing phase \(\phi\)
  • hbar (float) – (default 2) the value of \(\hbar\) in the commutation relation \([\x,\p]=i\hbar\)
Returns:

the squeezed state

Return type:

array

strawberryfields.utils.vacuum_state(basis='fock', fock_dim=5, hbar=2.0)[source]

Returns the vacuum state

Parameters:
  • basis (str) – If ‘fock’, calculates the initial state in the Fock basis. If ‘gaussian’, returns the vector of means and the covariance matrix.
  • fock_dim (int) – the size of the truncated Fock basis if using the Fock basis representation
  • hbar (float) – (default 2) the value of \(\hbar\) in the commutation relation \([\x,\p]=i\hbar\)
Returns:

the vacuum state

Return type:

array

strawberryfields.utils.coherent_state(a, basis='fock', fock_dim=5, hbar=2.0)[source]

Returns the coherent state

This can be returned either in the Fock basis,

\[|\alpha\rangle = e^{-|\alpha|^2/2} \sum_{n=0}^\infty \frac{\alpha^n}{\sqrt{n!}}|n\rangle\]

or as a Gaussian:

\[\mu = (\text{Re}(\alpha),\text{Im}(\alpha)),~~~\sigma = I\]

where \(\alpha\) is the displacement.

Parameters:
  • a (complex) – the displacement
  • basis (str) – If ‘fock’, calculates the initial state in the Fock basis. If ‘gaussian’, returns the vector of means and the covariance matrix.
  • fock_dim (int) – the size of the truncated Fock basis if using the Fock basis representation
  • hbar (float) – (default 2) the value of \(\hbar\) in the commutation relation \([\x,\p]=i\hbar\)
Returns:

the coherent state

Return type:

array

strawberryfields.utils.squeezed_state(r, p, basis='fock', fock_dim=5, hbar=2.0)[source]

Returns the squeezed state

This can be returned either in the Fock basis,

\[|z\rangle = \frac{1}{\sqrt{\cosh(r)}}\sum_{n=0}^\infty \frac{\sqrt{(2n)!}}{2^n n!}(-e^{i\phi}\tanh(r))^n|2n\rangle\]

or as a Gaussian:

\[\mu = (0,0)\]
\begin{align*} \sigma = R(\phi/2)\begin{bmatrix}e^{-2r} & 0 \\0 & e^{2r} \\\end{bmatrix}R(\phi/2)^T \end{align*}

where \(z = re^{i\phi}\) is the squeezing factor.

Parameters:
  • r (complex) – the squeezing magnitude
  • p (float) – the squeezing phase \(\phi\)
  • basis (str) – If ‘fock’, calculates the initial state in the Fock basis. If ‘gaussian’, returns the vector of means and the covariance matrix.
  • fock_dim (int) – the size of the truncated Fock basis if using the Fock basis representation
  • hbar (float) – (default 2) the value of \(\hbar\) in the commutation relation \([\x,\p]=i\hbar\)
Returns:

the squeezed state

Return type:

array

strawberryfields.utils.displaced_squeezed_state(a, r, phi, basis='fock', fock_dim=5, hbar=2.0)[source]

Returns the squeezed coherent state

This can be returned either in the Fock basis,

\[|\alpha,z\rangle = e^{-\frac{1}{2}|\alpha|^2-\frac{1}{2}{\alpha^*}^2 e^{i\phi}\tanh{(r)}} \sum_{n=0}^\infty\frac{\left[\frac{1}{2}e^{i\phi}\tanh(r)\right]^{n/2}}{\sqrt{n!\cosh(r)}} H_n\left[ \frac{\alpha\cosh(r)+\alpha^*e^{i\phi}\sinh(r)}{\sqrt{e^{i\phi}\sinh(2r)}} \right]|n\rangle\]

where \(H_n(x)\) is the Hermite polynomial, or as a Gaussian:

\[\mu = (\text{Re}(\alpha),\text{Im}(\alpha))\]
\begin{align*} \sigma = R(\phi/2)\begin{bmatrix}e^{-2r} & 0 \\0 & e^{2r} \\\end{bmatrix}R(\phi/2)^T \end{align*}

where \(z = re^{i\phi}\) is the squeezing factor and \(\alpha\) is the displacement.

Parameters:
  • a (complex) – the displacement
  • r (complex) – the squeezing magnitude
  • phi (float) – the squeezing phase \(\phi\)
  • basis (str) – If ‘fock’, calculates the initial state in the Fock basis. If ‘gaussian’, returns the vector of means and the covariance matrix.
  • fock_dim (int) – the size of the truncated Fock basis if using the Fock basis representation
  • hbar (float) – (default 2) the value of \(\hbar\) in the commutation relation \([\x,\p]=i\hbar\)
Returns:

the squeezed coherent state

Return type:

array

strawberryfields.utils.fock_state(n, fock_dim=5)[source]

Returns the Fock state

Parameters:
  • n (int) – the occupation number
  • fock_dim (int) – the size of the truncated Fock basis
Returns:

the Fock state

Return type:

array

strawberryfields.utils.cat_state(a, p=0, fock_dim=5)[source]

Returns the cat state

\[|cat\rangle = \frac{1}{\sqrt{2(1+e^{-2|\alpha|^2}\cos(\phi))}} \left(|\alpha\rangle +e^{i\phi}|-\alpha\rangle\right)\]

with the even cat state given for \(\phi=0\), and the odd cat state given for \(\phi=\pi\).

Parameters:
  • a (complex) – the displacement
  • p (float) – parity, where \(\phi=p\pi\). p=0 corresponds to an even cat state, and p=1 an odd cat state
  • fock_dim (int) – the size of the truncated Fock basis
Returns:

the cat state

Return type:

array

strawberryfields.utils.randnc(*arg)[source]

Normally distributed array of random complex numbers.

strawberryfields.utils.random_covariance(N, hbar=2, pure=False, block_diag=False)[source]

Random covariance matrix.

Parameters:
  • N (int) – number of modes
  • hbar (float) – the value of \(\hbar\) to use in the definition of the quadrature operators \(\x\) and \(\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 \(q\) do not mix with the momenta \(p\) and thus the covariance matrix is block diagonal.
Returns:

random \(2N\times 2N\) covariance matrix

Return type:

array

strawberryfields.utils.random_symplectic(N, passive=False, block_diag=False)[source]

Random symplectic matrix representing a Gaussian transformation.

The squeezing parameters \(r\) for active transformations are randomly sampled from the standard normal distribution, while passive transformations are randomly sampled from the Haar measure.

Parameters:
  • 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 \(q\) do not mix with the momenta \(p\) and thus the symplectic operator is block diagonal
Returns:

random \(2N\times 2N\) symplectic matrix

Return type:

array

strawberryfields.utils.random_interferometer(N, real=False)[source]

Random unitary matrix representing an interferometer.

For more details, see [43].

Parameters:
  • N (int) – number of modes
  • real (bool) – return a random real orthogonal matrix
Returns:

random \(N\times N\) unitary distributed with the Haar measure

Return type:

array

class strawberryfields.utils.operation(ns)[source]

Groups a sequence of gates into a single operation to be used within a Program context.

For example:

@sf.operation(3)
def custom_operation(v1, v2, q):
    CZgate(v1) | (q[0], q[1])
    Vgate(v2) | q[2]

Here, the operation decorator must recieve an argument detailing the number of subsystems the resulting custom operation acts on.

The function it acts on can contain arbitrary Python and Blackbird code that may normally be placed within a Program context. Note that it must always accept the register q it acts on as the last argument of the function.

Once defined, it can be used like any other quantum operation:

prog = sf.Program(3)
with prog.context as q:
    custom_operation(0.5719, 2.0603) | (q[0], q[1], q[3])

Note that here, we do not pass the register q directly to the function - instead, it is defined on the right hand side of the | operation, like all other Blackbird code.

Parameters:ns (int) – number of subsystems required by the operation
strawberryfields.utils.is_unitary(prog)[source]

True iff all the operations in the program are unitary.

Parameters:prog (Program) – quantum program
Returns:True iff all operations in the program are of type strawberryfields.ops.Gate
Return type:bool
strawberryfields.utils.is_channel(prog)[source]

True iff all the operations in the program can be represented as quantum channels.

Parameters:prog (Program) – quantum program
Returns:True if all operations in the program are of types strawberryfields.ops.Gate and strawberryfields.ops.Channel
Return type:bool
strawberryfields.utils.extract_unitary(prog, cutoff_dim: int, vectorize_modes: bool = False, backend: str = 'fock')[source]

Numerical array representation of a unitary quantum circuit.

Note that the circuit must only include operations of the strawberryfields.ops.Gate class.

  • If vectorize_modes=True, it returns a matrix.
  • If vectorize_modes=False, it returns an operator with \(2N\) indices, where N is the number of modes that the Program is created with. Adjacent indices correspond to output-input pairs of the same mode.

Example

This shows the Hong-Ou-Mandel effect by extracting the unitary of a 50/50 beamsplitter, and then computing the output given by one photon at each input (notice the order of the indices: \([out_1, in_1, out_2, in_2,\dots]\)). The result tells us that the two photons always emerge together from a random output port and never one per port.

>>> prog = sf.Program(num_subsystems=2)
>>> with prog.context as q:
>>>     BSgate(np.pi/4) | q
>>> U = extract_unitary(prog, cutoff_dim=3)
>>> print(abs(U[:,1,:,1])**2)
[[0.  0.  0.5]
 [0.  0.  0. ]
 [0.5 0.  0. ]])
Parameters:
  • prog (Program) – quantum program
  • cutoff_dim (int) – dimension of each index
  • vectorize_modes (bool) – If True, reshape input and output modes in order to return a matrix.
  • backend (str) – the backend to build the unitary; 'fock' (default) and 'tf' are supported
Returns:

numerical array of the unitary circuit

as a NumPy ndarray ('fock' backend) or as a TensorFlow Tensor ('tf' backend)

Return type:

array, tf.Tensor

Raises:

TypeError – if the operations used to construct the circuit are not all unitary

strawberryfields.utils.extract_channel(prog, cutoff_dim: int, representation: str = 'choi', vectorize_modes: bool = False)[source]

Numerical array representation of the channel corresponding to a quantum circuit.

The representation choices include the Choi state representation, the Liouville representation, and the Kraus representation.

Note

Channel extraction can currently only be performed using the 'fock' backend.

Tensor shapes

  • If vectorize_modes=True:
    • representation='choi' and representation='liouville' return an array with 4 indices
    • representation='kraus' returns an array of Kraus operators in matrix form
  • If vectorize_modes=False:
    • representation='choi' and representation='liouville' return an array with \(4N\) indices
    • representation='kraus' returns an array of Kraus operators with \(2N\) indices each, where \(N\) is the number of modes that the Program is created with

Note that the Kraus representation automatically returns only the non-zero Kraus operators. One can reduce the number of operators by discarding Kraus operators with small norm (thus approximating the channel).

Choi representation

Mathematically, the Choi representation of a channel is a bipartite state \(\Lambda_{AB}\) which contains a complete description of the channel. The way we use it to compute the action of the channel \(\mathcal{C}\) on an input state \(\mathcal{\rho}\) is as follows:

\[\mathcal{C}(\rho) = \mathrm{Tr}[(\rho_A^T\otimes\mathbb{1}_B)\Lambda_{AB}]\]

The indices of the non-vectorized Choi operator match exactly those of the state, so that the action of the channel can be computed as (e.g., for one mode or for vectorize_modes=True):

>>> rho_out = np.einsum('ab,abcd', rho_in, choi)

Notice that this respects the transpose operation.

For two modes:

>>> rho_out = np.einsum('abcd,abcdefgh', rho_in, choi)

Combining consecutive channels (in the order \(1,2,3,\dots\)) is also straightforward with the Choi operator:

>>> choi_combined = np.einsum('abcd,cdef,efgh', choi_1, choi_2, choi_3)

Liouville operator

The Liouville operator is a partial transpose of the Choi operator, such that the first half of consecutive index pairs are the output-input right modes (i.e., acting on the “bra” part of the state) and the second half are the output-input left modes (i.e., acting on the “ket” part of the state).

Therefore, the action of the Liouville operator (e.g., for one mode or for vectorize_modes=True) is

\[\mathcal{C}(\rho) = \mathrm{unvec}[\mathcal{L}\mathrm{vec}(\rho)]\]

where vec() and unvec() are the operations that stack the columns of a matrix to form a vector and vice versa. In code:

>>> rho_out = np.einsum('abcd,bd->ca', liouville, rho_in)

Notice that the state contracts with the second index of each pair and that we output the ket on the left (c) and the bra on the right (a).

For two modes we have:

>>> rho_out = np.einsum('abcdefgh,fbhd->eagc', liouville, rho_in)

The Liouville representation has the property that if the channel is unitary, the operator is separable. On the other hand, even if the channel were the identity, the Choi operator would correspond to a maximally entangled state.

The choi and liouville operators in matrix form (i.e., with two indices) can be found as follows, where D is the dimension of each vectorized index (i.e., for \(N\) modes, D=cutoff_dim**N):

>>> choi_matrix = liouville.reshape(D**2, D**2).T
>>> liouville_matrix = choi.reshape(D**2, D**2).T

Kraus representation

The Kraus representation is perhaps the most well known:

\[\mathcal{C}(\rho) = \sum_k A_k\rho A_k^\dagger\]

So to define a channel in the Kraus representation one needs to supply a list of Kraus operators \(\{A_k\}\). In fact, the result of extract_channel in the Kraus representation is a rank-3 tensor, where the first index is the one indexing the list of operators.

Adjacent indices of each Kraus operator correspond to output-input pairs of the same mode, so the action of the channel can be written as (here for one mode or for vectorize_modes=True):

>>> rho_out = np.einsum('abc,cd,aed->be', kraus, rho_in, np.conj(kraus))

Notice the transpose on the third index string (aed rather than ade), as the last operator should be the conjugate transpose of the first, and we cannot just do np.conj(kraus).T because kraus has 3 indices and we just need to transpose the last two.

Example

Here we show that the Choi operator of the identity channel is proportional to a maximally entangled Bell \(\ket{\phi^+}\) state:

>>> prog = sf.Program(num_subsystems=1)
>>> C = extract_channel(prog, cutoff_dim=2, representation='choi')
>>> print(abs(C).reshape((4,4)))
[[1. 0. 0. 1.]
 [0. 0. 0. 0.]
 [0. 0. 0. 0.]
 [1. 0. 0. 1.]]
Parameters:
  • prog (Program) – program containing the circuit
  • cutoff_dim (int) – dimension of each index
  • representation (str) – choice between 'choi', 'liouville' or 'kraus'
  • vectorize_modes (bool) – if True, reshapes the result into rank-4 tensor, otherwise it returns a rank-4N tensor, where N is the number of modes
Returns:

channel, according to the specified options

Return type:

array

Raises:

TypeError – if the gates used to construct the circuit are not all unitary or channels