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 mode value.
mag(x) Returns the magnitude \(|z|\) of a measured mode value.
phase(x) Returns the phase \(\phi\) of a measured mode value \(z=re^{i\phi}\).
scale(x, a) Scales the value of a measured mode by factor a.
shift(x, b) Shifts the value of a measured mode by factor b.
scale_shift(x, a, b) Scales the value of a measured mode by factor a then shifts the result by b.
power(x, a) Raises the value of a measured mode 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]) Returns a random covariance matrix.
random_symplectic(N[, passive]) Returns a random symplectic matrix representing a Gaussian transformation.
random_interferometer(N) Returns a 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 an engine context.

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

Engine functions

These functions act on instances of a compiler Engine, returning or extracting information from the queued engine operations.

For example, these might be used directly on the engine object as follows:

eng, q = sf.Engine(2)

with eng:
    BSgate(0.543, 0.123) | (q[0], q[1])

U = extract_unitary(eng, cutoff_dim=10)

In this example, U is an array representing the unitary operation in the Fock basis of the queued engine operations (here, a single beamsplitter).

is_unitary(engine) Returns True if every command in the queue is of type strawberryfields.ops.Gate.
is_channel(engine) Returns true if every command in the queue is either of type strawberryfields.ops.Gate or type strawberryfields.ops.Channel.
extract_unitary(engine, cutoff_dim, …) Returns the array representation of a queued unitary circuit as a NumPy ndarray ('fock' backend) or as a TensorFlow Tensor ('tf' backend).
extract_channel(engine, cutoff_dim, …) Returns a numerical array representation of a channel from a queued Engine circuit.

Note

These functions act on an engine with queued operations. If the engine has already been run via eng.run(), then the queue will be empty and these functions will simply be acting on an empty circuit (i.e., the identity).

Code details

strawberryfields.utils.neg(x)[source]

Negates a measured mode value.

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

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

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

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

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

Scales the value of a measured mode by factor a.

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

Shifts the value of a measured mode 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 value of a measured mode 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 value of a measured mode 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)[source]

Returns a 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
Returns:

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

Return type:

array

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

Returns a 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.
Returns:

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

Return type:

array

strawberryfields.utils.random_interferometer(N)[source]

Returns a random unitary matrix representing an interferometer.

For more details, see [41].

Parameters:N (int) – number of modes
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 an engine 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 an engine context. Note that it must always accept the qumode register q it acts on as the last argument of the function.

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

eng, q = sf.Engine(3)
with eng:
    custom_operation(0.5719, 2.0603) | (q[0], q[1], q[3])

Note that here, we do not pass the qumode 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 registers required by the operation
strawberryfields.utils.is_unitary(engine)[source]

Returns True if every command in the queue is of type strawberryfields.ops.Gate.

Parameters:engine (Engine) – a Strawberry Fields engine instance
Returns:returns True if all queued operations are unitary
Return type:bool
strawberryfields.utils.is_channel(engine)[source]

Returns true if every command in the queue is either of type strawberryfields.ops.Gate or type strawberryfields.ops.Channel.

Parameters:engine (Engine) – a Strawberry Fields engine instance
Returns:returns True if all queued operations are either a quantum gate or a channel
Return type:bool
strawberryfields.utils.extract_unitary(engine, cutoff_dim: int, vectorize_modes: bool = False, backend: str = 'fock')[source]

Returns the array representation of a queued unitary circuit as a NumPy ndarray ('fock' backend) or as a TensorFlow Tensor ('tf' backend).

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 engine 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.

>>> engine, (A, B) = sf.Engine(num_subsystems=2)
>>> with engine:
>>>     BSgate(np.pi/4) | (A, B)
>>> U = extract_unitary(engine, cutoff_dim=3)
>>> print(abs(U[:,1,:,1])**2)
[[0.  0.  0.5]
 [0.  0.  0. ]
 [0.5 0.  0. ]])
Parameters:
  • engine (Engine) – the engine containing a queued circuit
  • 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:

the numerical array of the unitary circuit

Return type:

array

Raises:

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

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

Returns a numerical array representation of a channel from a queued Engine 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 engine 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:

>>> engine, A = sf.Engine(num_subsystems=1)
>>> C = extract_channel(engine, 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:
  • engine (Engine) – the engine 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:

the numerical array of the channel, according to the specified representation and vectorization options

Raises:

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