Getting started

Background

The Walrus API

  1. Docs
  2. The Walrus 0.22.0-dev documentation

The Walrus Documentation

Release:

0.22.0-dev

Responsive image

A library for the calculation of hafnians, Hermite polynomials, and Gaussian boson sampling.

Using The Walrus

See the quick guide for an overview of available functions in The Walrus

Background

Learn about the hafnian, loop hafnian, and its relationship to quantum photonics

API

Explore The Walrus Python API

Features

  • Fast calculation of hafnians, loop hafnians, and torontonians of general and certain structured matrices powered by Numba.

  • An easy to use interface to use the loop hafnian for Gaussian quantum state calculations.

  • Sampling algorithms for hafnian and torontonians of graphs.

  • Efficient classical methods for approximating the hafnian of non-negative matrices.

  • Easy to use implementations of the multidimensional Hermite polynomials, which can also be used to calculate hafnians of all reductions of a given matrix.

Getting started

To get the The Walrus installed and running on your system, begin at the download and installation guide. Then, familiarise yourself with some background information on the Hafnian and the computational algorithm.

For getting started with using the The Walrus in your own code, have a look at the Python tutorial.

Finally, detailed documentation on the code and API is provided.

Support

If you are having issues, please let us know, either by email or by posting the issue on our Github issue tracker.

Authors

Nicolás Quesada, Brajesh Gupt, and Josh Izaac.

If you are doing research using The Walrus, please cite our paper:

Brajesh Gupt, Josh Izaac and Nicolás Quesada. The Walrus: a library for the calculation of hafnians, Hermite polynomials and Gaussian boson sampling. Journal of Open Source Software, 4(44), 1705 (2019)

License

The Walrus library is free and open source, released under the Apache License, Version 2.0.

Installation and Downloads

The Walrus requires Python version 3.7, 3.8, 3.9, or 3.10. Installation of The Walrus, as well as all dependencies, can be done using pip:

pip install thewalrus

Compiling from source

The Walrus has the following dependencies:

You can compile the latest development version by cloning the git repository, and installing using pip in development mode.

$ git clone https://github.com/XanaduAI/thewalrus.git
$ cd thewalrus && python -m pip install -e .

Software tests

To ensure that The Walrus library is working correctly after installation, the test suite can be run locally using pytest.

Additional packages are required to run the tests. These dependencies can be found in requirements-dev.txt and can be installed using pip:

$ pip install -r requirements-dev.txt

To run the tests, navigate to the source code folder and run the command

$ make test

Documentation

The Walrus documentation is available online on Read the Docs.

Additional packages are required to build the documentation locally as specified in doc/requirements.txt. These packages can be installed using:

$ sudo apt install pandoc
$ pip install -r docs/requirements.txt

To build the HTML documentation, go to the top-level directory and run the command

$ make doc

The documentation can then be found in the docs/_build/html/ directory.

Research and contribution

Research

If you are doing research using The Walrus, please cite our paper:

Brajesh Gupt, Josh Izaac and Nicolás Quesada. The Walrus: a library for the calculation of hafnians, Hermite polynomials and Gaussian boson sampling. Journal of Open Source Software, 4(44), 1705 (2019)

We are always open for collaboration, and can be contacted at research@xanadu.ai.

Contribution

If you would like to get involved with The Walrus, or you would like to contribute, simply fork our Github repository, and make a descriptive pull request.

Support

If you are having issues, please let us know, either by email or by posting the issue on our Github issue tracker.

We have a mailing list located at: support@xanadu.ai

Quick guide

This section provides a quick guide to find which function does what in The Walrus.

What you want

Corresponding function

Numerical hafnian

thewalrus.hafnian()

Symbolic hafnian

thewalrus.reference.hafnian()

Hafnian of a matrix with repeated rows and columns

thewalrus.hafnian_repeated()

Hafnians of all possible reductions of a given matrix

thewalrus.hafnian_batched()

Hafnian samples of a Gaussian state

thewalrus.samples.hafnian_sample_state()

Torontonian samples of a Gaussian state

thewalrus.samples.torontonian_sample_state()

Hafnian samples of a graph

thewalrus.samples.hafnian_sample_graph()

Torontonian samples of a graph

thewalrus.samples.torontonian_sample_graph()

All probability amplitudes of a pure Gaussian state

thewalrus.quantum.state_vector()

All matrix elements of a general Gaussian state

thewalrus.quantum.density_matrix()

A particular probability amplitude of pure Gaussian state

thewalrus.quantum.pure_state_amplitude()

A particular matrix element of a general Gaussian state

thewalrus.quantum.density_matrix_element()

The Fock representation of a Gaussian unitary

thewalrus.quantum.fock_tensor()

Note that all the hafnian functions listed above generalize to loop hafnians.

The hafnian

Section author: Nicolás Quesada <nicolas@xanadu.ai>

The hafnian

The hafnian of an \(n \times n\), symmetric matrix \(\bm{A} =\bm{A}^T\) is defined as [1]

\[\label{eq:hafA} \haf(\bm{A}) = \sum_{M \in \text{PMP}(n)} \prod_{\scriptscriptstyle (i, j) \in M} A_{i, j},\]

where \(\text{PMP}(n)\) stands for the set of perfect matching permutations of \(n\) (even) objects, i.e., permutations \(\sigma:[n]\rightarrow [n]\) such that \(\forall i:\sigma(2i-1)<\sigma(2i)\) and \(\forall i:\sigma(2i-1)<\sigma(2i+1)\).

It was so named by Eduardo R. Caianiello [2] “to mark the fruitful period of stay in Copenhagen (Hafnia in Latin).” [3]

The set PMP(\(n\)) used to define the hafnian contains

\[\label{eq:haf1} |\text{PMP}(n)|=(n-1)!! = 1 \times 3 \times 5 \times \ldots \times (n -1),\]

elements and thus as defined it takes \((n-1)!!\) additions of products of \(n/2\) numbers to calculate the hafnian of an \(n \times n\) matrix. Note that the diagonal elements of the matrix \(\bm{A}\) do not appear in the calculation of the hafnian and are (conventionally) taken to be zero.

For \(n=4\) the set of perfect matchings is

\[\label{eq:PMP4} \text{PMP}(4) = \big\{ (0,1)(2,3),\ (0,2)(1,3),\ (0,3),(1,2) \big\},\]

and the hafnian of a \(4 \times 4\) matrix \(\bm{B}\) is

\[\label{eq:hafB} \haf(\bm{B}) = B_{0,1} B_{2,3}+B_{0,2}B_{1,3}+B_{0,3} B_{1,2}.\]

The hafnian of an odd sized matrix is defined to be zero; if \(\bm{A}=\bm{A}^T\) is \(M\) dimensional and \(M\) is odd then \(\haf(\bm{A}) = 0\). Note that, for convenience, we define the hafnian of an empty matrix, i.e., a matrix of dimension zero by zero, to be 1.

The hafnian is a homogeneous function of degree \(n/2\) in the matrix entries of an \(n \times n\) matrix \(\bm{A}\). This implies that

\[\haf(\mu \bm{A}) = \mu ^{n/2} \haf(\bm{A}),\]

where \(\mu\) is a scalar. More generally, if \(\bm{W} = \text{diag}(w_0,\ldots,w_{n-1})\), then it holds that (see proposition 4.2.3 of [1])

\[\haf( \bm{W} \bm{A} \bm{W} ) = \left(\prod_{i=0}^{n-1} w_i\right) \haf(\bm{A}).\]

The definition used to introduce the hafnian is rather algebraic and brings little intuition. To gain more insight in the next section we introduce some graph theory language and use it to present a more intuitive vision of what the hafnian “counts”.

Basics of graphs

A graph is an ordered pair \(G=(V,E)\) where \(V\) is the set of vertices, and \(E\) is the set of edges linking the vertices of the graph, i.e., if \(e \in E\) then \(e=(i,j)\) where \(i,j \in V\). In this section we consider graphs without loops (we relax this in the next section), thus we do not allow for edges \(e = (i,i)\) connecting a given vertex to itself. A 6 vertex graph is shown here

_images/graph.svg

the vertices are labelled \(V = \{1,2,3,4,5,6 \}\) and the edges are \(E=\{(1,4),(2,4),(2,5),(3,4),(3,5),(3,6),(5,5) \}\).

A matching \(M\) is a subset of the edges in which no two edges share a vertex. An example of matching is \(M=(1,4)(3,6)\) represented by the blue lines in the following figure

_images/matching.svg

In the figure above we know we have a matching because none of the highlighted edges shares a vertex.

A perfect matching is a matching which matches all the vertices of the graph, such as for example \(M=(1,4)(2,5)(3,6)\), which is represented again by the blue lines in the following figure

_images/pm.svg

The blue lines represent a perfect matching because, they are a matching, i.e., the edges do no overlap on any vertex and all the vertices are covered by one and only edge.

A complete graph is a graph where every vertex is connected to every other vertex. For loopless graphs having \(n\) vertices, the number of perfect matchings is precisely [1]

\[|\text{PMP}(n)|=(n-1)!! = 1 \times 3 \times \ldots \times (n-1).\]

where we use \(\text{PMP}(n)\) to indicate the set of perfect matchings of introduced in the previous section, and the notation \(|V|\) to indicate the number of elements in the set \(V\). Note that this number is nonzero only for even \(n\), since for odd \(n\) there will always be one unmatched vertex.

In the following figure we illustrate the 3 perfect matchings of a complete graph with 4 vertices

_images/pmp4.svg

Perfect matchings and hafnians

An important question concerning a given graph \(G=(V,E)\) is the number of perfect matchings it has. One possible way to answer this question is to iterate over the perfect matchings of a complete graph and at each step check if the given perfect matching of the complete graph is also a perfect matching of the given graph. A simple way to automatize this process is by constructing the adjacency matrix of the graph. The adjacency matrix \(\bm{A}\) of a graph \(G=(V,E)\) is a 0-1 matrix that has \(\bm{A}_{i,j} = \bm{A}_{j,i}=1\) if, and only if, \((i,j) \in E\) and 0 otherwise. For the example graph in the previous section, the adjacency matrix is

\[\begin{split}\bm{A}' = \begin{bmatrix} 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 1 & 0 \\ 0 & 0 & 0 & 1 & 1 & 1 \\ 1 & 1 & 1 & 0 & 0 & 0 \\ 0 & 1 & 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 \end{bmatrix}.\end{split}\]

The number of perfect matchings of a (loopless) graph is simply given by the hafnian of its adjacency matrix

\[\text{haf}(\bm{A}) = \sum_{M \in \text{PMP}(n)} \prod_{\scriptscriptstyle (i,j) \in M} {A}_{i,j}.\]

For the graph in the previous section we can easily confirm that the perfect matching we found is the only perfect matching since

\[\text{haf}(\bm{A}') = 1.\]

The definition of the hafnian immediately generalizes to weighted graphs, where we assign a real or complex number to the entries of the symmetric matrix \(\bm{A}\).

Special values of the hafnian

Here we list some special values of the hafnian for certain special matrices.

  • If the matrix \(\bm{A}\) has the following block form

\[\begin{split}\bm{A}_{\text{block}} = \left[\begin{array}{c|c} 0 & \bm{C} \\ \hline \bm{C}^T & 0 \\ \end{array}\right],\end{split}\]

then it holds that \(\text{haf}\left( \bm{A}_{\text{block}} \right) = \text{per}(\bm{C})\) where \(\text{per}\) is the permanent matrix function defined as [1]

\[\text{per}(\bm{C})=\sum_{\sigma\in S_n}\prod_{i=1}^n C_{i,\sigma(i)}.\]

The sum here extends over all elements \(\sigma\) of the symmetric group \(S_n\).

  • If \(\bm{A}_{\text{rank-one}} = \bm{e} \bm{e}^T\) is a rank one matrix of size \(n\) then

\[\text{haf}\left( \bm{A}_{\text{rank-one}} \right) = (n-1)!! \prod_{i=1}^{n-1} e_i.\]

In particular, the hafnian of the all ones matrix is precisely \((n-1)!!\).

  • If \(\bm{A}_{\text{direct sum}} = \bm{A}_1 \oplus \bm{A}_2\) is a block diagonal matrix then

\[\text{haf}\left(\bm{A}_{\text{direct sum}}\right) = \text{haf}\left( \bm{A}_1 \oplus \bm{A}_2 \right) = \text{haf}\left( \bm{A}_1 \right) \text{haf}\left( \bm{A}_2 \right).\]

This identity simply expresses the fact that the number of perfect matchings of a graph that is made of two disjoint subgraphs is simply the product of the number of perfect matchings of the two disjoint subgraphs.

The loop hafnian

Section author: Nicolás Quesada <nicolas@xanadu.ai>

Graphs with loops

In the previous section we introduced the hafnian as a way of counting the number of perfect matchings of a loopless graph. The loop hafnian does the same for graphs with loops. Before defining the loop hafnian let us introduce graphs with loops. A graph will still be an ordered pair \((V,E)\) of vertices and edges but now we will allow edges of the form \((i,i)\) where \(i \in V\). We can now define adjacency matrices in the same way as we did before i.e. if \((i,j) \in E\) then \(M_{i,j}=1\) and otherwise \(M_{i,j}=0\).

Consider the following graph

_images/loop_fig1.svg

for which we have \(V = \{1,2,3,4,5,6 \}\), the edges are \(E=\{(1,1),(1,4),(2,4),(2,5),(3,4),(3,5),(3,6),(5,5) \}\) and the adjacency matrix is

\[\begin{split}\bm{A}'' = \begin{bmatrix} 1 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 1 & 0 \\ 0 & 0 & 0 & 1 & 1 & 1 \\ 1 & 1 & 1 & 0 & 0 & 0 \\ 0 & 1 & 1 & 0 & 1 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 \end{bmatrix}.\end{split}\]

Note that there are now nonzero elements in the diagonal indicating that vertices 1 and 5 have loops.

Once we allow for loops we have more options for making perfect matchings. For example for the graph shown above there are now 2 perfect matchings, illustrated in blue in the following figure

_images/loop_fig2.svg

As was done before for the hafnian we introduce the set of single pair matchings \(\text{SPM}(n)\) as the set of perfect matchings of a graph of size \(n`with loops :cite:`bjorklund2018faster\). For \(n=4\) we have

\[\begin{split}\text{SPM}(4) = \big\{ (0,1)(2,3),\ (0,2)(1,3),\ (0,3),(1,2),\ (0,0)(1,1)(2,3), \ (0,1)(2,2)(3,3),\\ (0,0)(2,2)(1,3),\ (0,2)(1,1)(3,3),\ (0,0)(3,3)(1,2),\ (0,3)(1,1)(2,2),\ (0,0)(1,1)(2,2)(3,3)\big\}.\end{split}\]

For a graph with 4 vertices they are

_images/fig2.svg

Note that there is a one to one correspondence (a bijection) between the elements in \(\text{SPM}(n)\) and the number of matchings of a graph with \(n\) vertices \(H(n)\). To see why this is the case, note that any element of \(\text{SPM}(n)\) can be converted into a matching by removing all the vertices that are loops. For example, to the following element \((0,0)(2,2)(1,3)\) we associate the matching \((1,2)\). Note that this mapping is one-to-one since, given a matching, we can always add as loops all the other vertices that are not part of the matching. Using this bijection we conclude that the number of elements in \(\text{SPM}(n)\) is (see The On-Line Encyclopedia of Integer Sequences)

\[|\text{SPM}(n)| = T(n),\]

where \(T(n)\) is the \(n^{\text{th}}\) telephone number.

Note that in general for given graph size \(n\) there a lot more single pair matching that there are perfect matchings. Their ratio goes like [4]

\[\frac{\text{SPM}(n)}{\text{PMP}(n)} = \frac{T(n)}{(n-1)!!} \sim e^{\sqrt{n}}.\]

The loop hafnian

We will also be interested in a generalization of the hafnian function where we will now allow for adjacency matrices that have loops. This new function we call the loop hafnian (lhaf). As explained before, the weight associated with said loops will be allocated in the diagonal elements of the adjacency matrix \(\bm{A}\) (which were previously ignored in the definition of the hafnian). To account for the possibility of loops we generalized the set of perfect matching permutations PMP to the single-pair matchings (SPM). Thus we define [4]

\[\lhaf(\bm{A}) = \sum_{M \in \text{SPM}(n)} \prod_{\scriptscriptstyle (i,j) \in M} A_{i,j}.\]

The lhaf of a \(4 \times 4\) matrix \(\bm{B}\) is

\[\begin{split}\lhaf(\bm{B}) =& B_{0,1} B_{2,3}+B_{0,2}B_{1,3}+B_{0,3} B_{1,2}\\ &+ B_{0,0} B_{1,1} B_{2,3}+B_{0,1} B_{2,2} B_{3,3}+B_{0,2}B_{1,1}B_{3,3}\nonumber\\ &+ B_{0,0} B_{2,2} B_{1,3}+B_{0,0}B_{3,3}B_{1,2}+B_{0,3} B_{1,1} B_{2,2}\nonumber\\ &+ B_{0,0} B_{1,1} B_{2,2} B_{3,3}. \nonumber\end{split}\]

Finally, let us comment on the scaling properties of the \(\haf\) and \(\lhaf\). Unlike the hafnian, the loop hafnian is not homogeneous in its matrix entries, i.e.

\[\begin{split}\haf(\mu \bm{A}) &= \mu ^{n/2} \haf(\bm{A}) \text{ but},\\ \lhaf(\mu \bm{A}) &\neq \mu ^{n/2} \lhaf(\bm{A}).\end{split}\]

where \(n\) is the size of the matrix \(\bm{A}\) and \(\mu \geq 0\). However if we split the matrix \(\bm{A}\) in terms of its diagonal \(\bm{A}_{\text{diag}}\) part and its offdiagonal part \(\bm{A}_{\text{off-diag}}\)

\[\bm{A} = \bm{A}_{\text{diag}}+\bm{A}_{\text{off-diag}},\]

then it holds that [4]

\[\lhaf(\sqrt{\mu} \bm{A}_{\text{diag}}+ \mu \bm{A}_{\text{off-diag}}) = \mu^{n/2} \lhaf(\bm{A}_{\text{diag}}+ \bm{A}_{\text{off-diag}}) =\mu^{n/2} \lhaf(\bm{A}).\]

One can use the loop hafnian to count the number of matchings of a loopless graph by simply calculating the loop hafnian of its adjacency matrix adding ones in its diagonal.

Finally, if \(\bm{A}_{\text{direct sum}} = \bm{A}_1 \oplus \bm{A}_2\) is a block diagonal matrix then

\[\text{lhaf}\left(\bm{A}_{\text{direct sum}}\right) = \text{lhaf}\left( \bm{A}_1 \oplus \bm{A}_2 \right) = \text{lhaf}\left( \bm{A}_1 \right) \text{lhaf}\left( \bm{A}_2 \right).\]

As for the hafnian, this identity tell us that the number of perfect matchings of a graph that is made of two disjoint subgraphs is simply the product of the number of perfect matchings of the two disjoint subgraphs.

The algorithms

Section author: Nicolás Quesada <nicolas@xanadu.ai>

The exact calculation of the number of perfect matchings (and the hafnian) for general graphs (symmetric matrices) has been investigated by several authors in recent years. An algorithm running in \(O(n^2 2^n)\) time was given by Björklund and Husfeldt [5] in 2008. In the same paper an algorithm running in \(O(1.733^n)\) time was presented using fast matrix multiplication. In the same year Kan [6] presented an algorithm running in time \(O(n 2^n)\) by representing the hafnian as a moment of the multinormal distribution. Koivisto [7] gave an \(O^*(\phi^n)\) time and space algorithm, where \(\phi = (1+\sqrt{5})/2 \approx 1.618\) is the Golden ratio. Nederlof [8] provided a polynomial space algorithm running in \(O(1.942^n)\) time.

Finally, Björklund [9] [4] and later Cygan and Pilipczuk [10] provided \(O(\text{poly}(n) 2^{n/2})\) time and polynomial space algorithms for the calculation of the general ring hafnian. These algorithms are believed to be close to optimal unless there are surprisingly efficient algorithms for the Permanent. This is because these two algorithms can also be used to count (up to polynomial corrections) the number of perfect matchings for bipartite graphs with the same exponential growth as Ryser’s algorithm for the permanent [11]. Equivalently, if one could construct an algorithm that calculates hafnians in time \(O(\alpha^{n/2})\) with \(\alpha<2\) one could calculate permanents faster than Ryser’s algorithm (which is the fastest known algorithm to calculate the permanent [12]). This is because of the identity

\[\begin{split}\text{haf} \left( \left[ \begin{array}{cc} 0 & \bm{W} \\ \bm{W}^T & 0 \\ \end{array} \right]\right) = \text{per}(\bm{W}),\end{split}\]

which states that a bipartite graph with two parts having \(n/2\) elements can always be thought as a simple graph with \(n\) vertices. It should be noted that improving over Ryser’s algorithm is a well-known open problem: e.g. Knuth [13] asks for an arithmetic circuit for the permanent with less than \(2^n\) operations. Also note that since the exact calculation of the permanent of \((0,1)\) matrices is in the #P complete class [14] the above identity shows that deciding if the hafnian of a complex matrix is larger than a given value is also in the #P complete class.

Tip

Permanents can be calculated directly using Ryser’s algorithm via thewalrus.perm().

Finally, note that the approximate methods developed for counting perfect matchings are aimed at (weighted-)graphs with real or positive entries [15, 16]. Of particular note is the approximate algorithm introduced by Barvinok for matrices with non-negative entries [17] further analyzed in Ref. [18].

In what follows we provide a pseudo-code or equations that give a basic intuition for the algorithms implemented in this library. The reader is referred to the original literature for proof of correctness and complexity.

Reference algorithm

We provide a reference implementation of the hafnian and loop hafnian that iterates over the sets \(\text{PMP}(n)\) and \(\text{SPM}(n)\). These implementations are extremely slow for even moderate sized matrices and are only provided for educational purposes.

Tip

Implemented as thewalrus.reference.hafnian(). The optional argument loops=True can be used to calculate loop hafnians.

Recursive algorithm

In 2012 Björklund [9] introduced the following algorithm to calculate the hafnian of a matrix of size \(n\) (even) in any field \(\mathbb{F}\) in time \(O(n^4 \log(n) 2^{n/2})\)

_images/bjorklund.svg

In the pseudocode above the following notation is used:

  • \([n]=\{0,1,2,\ldots,n-1\}\) is the set of the first \(n\) integers,

  • \(|X|\) is used to denote the number of elements in the set \(X\), and

  • \(P(X)\) is used to denote the power set, which is the set of all the subsets of the set \(X\). Note that if \(X\) has \(|X|\) elements, then its power set has \(2^{|X|}\) elements.

Note that the subindices and superindices in the matrices \(\bm{B}\) are not used for components of the matrices but rather to denote stages in the computation. Furthermore, these matrices contain polynomials in the symbolic variable \(r\) and that the final answer is obtained by adding the coefficients of \(r^{n/2}\) in the polynomial \(g\) at each step. In the implementation provided in this library the algorithm sketched above in pseudocode is turned into a recursion relation, hence the name we give it here.

Unfortunately, there is no known generalization of this algorithm to loop hafnians.

Tip

Implemented as thewalrus.hafnian(). This is the default algorithm for calculating hafnians.

Trace algorithm

Based on the work of Cygan and Pilipczuk [10], Björklund et al [4] introduced another algorithm to calculate the hafnian of a real or complex matrix of size \(n\) in 2018. This algorithm which runs in time \(O(n^3 2^{n/2})\) and can be more succinctly expressed as an equation

\[\text{haf}(\bm{A}) = \sum_{S \in P([n/2])} (-1)^{ |S|} f\left((\bm{A} \bm{X})_{S}\right),\]

where the matrix \(\bm{X}\) is defined as

\[\begin{split}\bm{X}= \bm{X}^T=\bm{X}^{-1}= \begin{bmatrix} \bm{0} & \mathbb{I} \\ \mathbb{I} & \bm{0} \end{bmatrix},\end{split}\]

\(\mathbb{I}\) is the identity matrix and the function \(f(\bm{C})\) takes a matrix \(\bm{C}\) and returns the coefficient of \(\eta^{n/2}\) in the following polynomial:

\[p_{n/2}(\eta \bm{C}) = \sum_{j=1}^{n/2} \frac{1}{j!}\left(\sum_{k=1}^{n/2} \frac{\text{tr}(\bm{C}^k)}{2k}\eta^k \right)^j.\]

This coefficient can be found by taking derivatives [19]

\[f(\bm{C}) = \frac{1}{(n/2)!} \left. \frac{d^{n/2}}{d\eta^{n/2}} p_{n/2}(\eta \bm{C}) \right|_{\eta=0} = \frac{1}{(n/2)!} \left. \frac{d^{n/2}}{d\eta^{n/2}} \frac{1}{\sqrt{\det(\mathbb{I} - \eta \bm{C})}} \right|_{\eta=0}.\]

The function \(p_{n/2}(\eta\bm{C})\) requires only the traces of the matrix powers of the matrices \(\bm{C}^k\), hence the name of this algorithm. These powertraces can be calculated using the characteristic polynomial of the input matrix using the La Budde algorithm [20].

This formula generalizes to the loop hafnian as follows

\[\text{lhaf}(\bm{A}) = \sum_{S \in P([n/2])} (-1)^{ |S|} f\left((\bm{A} \bm{X})_{S}\right),\]

where now the function \(f(\bm{C})\) takes a matrix \(\bm{C}\) and returns the coefficient of \(\eta^{n/2}\) in the following polynomial:

\[q(\eta, \bm{C}, \bm{v} ) = \sum_{j=1}^{n/2} \frac{1}{j!} \left(\sum_{k=1}^{n/2} \left( \frac{\text{Tr}(\bm{C}^k) }{(2k)} +\frac{\bm{v} (\bm{X} \bm{B})^{k-1} \bm{v}^T}{2} \right) \eta^k \right)^j.\]

where \(\bm{v} = \text{diag}(\bm{C})\) and we used the reduction operation (cf. notation) in terms of the set \(S\).

Tip

Implemented as thewalrus.hafnian() with the argument recursive=False. The loop hafnian calculation can be done by setting the option loops=True.

Repeated-moment algorithm

By mapping the calculation of moments of the multinormal distribution to the calculation of the hafnian, Kan [6] derived the following equation for the loop hafnian

\[\begin{split}\text{lhaf}\left( \text{vid}(\bm{B}_\bm{m},\bm{u}_\bm{m}) \right) &= \sum_{\nu_0=0}^{m_0} \ldots \sum_{\nu_{n-1}}^{m_{n-1}} \sum_{r=0}^{[m/2]} (-1)^{\sum_{i=0}^{n-1} \nu_i} {m_0 \choose \nu_0} \ldots {m_{n-1} \choose \nu_{n-1}} \frac{\left( \frac{\bm{h}^T \bm{B} \ \bm{h}}{2} \right)^r \left(\bm{h}^T \bm{u} \right)^{m-2r}}{r! (m-2r)!}, \\ \bm{h} &= \left(\tfrac{m_{0}}{2}-\nu_0,\ldots,\tfrac{m_{n-1}}{2}-\nu_{n-1} \right),\\ m&=m_0+\ldots+m_{n-1},\end{split}\]

where we used the reduction and vector in diagonal (\(\text{vid}\)) operations introduced in the notation section.

Note that if we pick \(m_i=1 \ \forall i\) and \(\bm{u} = \text{diag}(\bm{A})\) we recover the loop hafnian of \(\bm{A}\). In this case, the calculation of the loop hafnian requires \(O(n 2^n)\) operations, which is quadratically worse than Björklund’s algorithms. This formula is however useful when calculating hafnians and loop hafnians of matrices with repeated rows and columns for which column and row \(i\) are repeated \(m_i\) times; taking only \(O(n A G^n)\) operations to calculate the loop hafnian, where

\[\begin{split}A &= \frac{1}{n} \sum_{i=0}^{n-1} (m_i+1), \\ G &= \left( \prod_{i=0}^{n-1}(m_i+1) \right)^{1/n}.\end{split}\]

Compare this with Björklund’s algorithm, which requires \(O\left((A n)^3 \left(\sqrt{2}^{A}\right)^n\right)\) operations.

Tip

Implemented as thewalrus.hafnian_repeated(). The vector \(\bm{m}\) is passed in the variable rpt. The loop hafnian calculation can be done by passing the variable mu with the values of the vector \(\bm{u}\) and the option loops=True.

Batched algorithm

Using the multidimensional Hermite polynomials, and their connection to the matrix elements of Gaussian states and hafnians discussed in the next section, one can calculate the hafnians of all reductions of a matrix \(\bm{B} \in \mathbb{C}^{n \times n}\) up to a given cutoff. The reduction of matrix \(\bm{B}\) is precisely the matrix \(\bm{B}_{\bm{m}}\) obtained by repeating (or removing) the \(i^{\text{th}}\) row and column \(m_i\) times. Thus given a cutoff \(m_{\max}\), one can use the batched algorithm to calculate

\[\text{haf}\left( \bm{B}_\bm{m} \right)\]

for all \(0\leq m_i < m_{\max}\), thus this function returns a tensor with \((m_{\max})^n\) components.

One can also use this function to calculate the same loop hafnians that Kan’s algorithm returns

\[\text{lhaf}\left( \text{vid}(\bm{B}_\bm{m},\bm{u}_\bm{m}) \right)\]

if provided also with a vector \(\bm{u}\). Note that this parameter is optional.

Internally, these hafnians are calculated by using the recursion relation of the multidimensional Hermite polynomials discussed in the next section.

Tip

Implemented as thewalrus.hafnian_batched(). The loop hafnian calculation can be done by passing the variable mu with the values of the vector \(\bm{u}\).

Approximate algorithm

In 1999 Barvinok [17] provided a surprisingly simple algorithm to approximate the hafnian of a symmetric matrix with non-negative entries. Let the matrix have entries \(A_{i,j}\) and define the antisymmetric stochastic matrix with entries that distribute according to \(W_{i,j} = -W_{i,j} \sim \mathcal{N}(0,A_{i,j})\), where \(\mathcal{N}(\mu,\sigma^2)\) is the normal distribution with mean \(\mu\) and variance \(\sigma^2\). The following now holds:

\[\text{haf}(\bm{A}) = \mathbb{E} \left( \text{det}(\bm{W}) \right)\]

where \(\mathbb{E}\) denotes the usual statistical expectation value, and \(\text{det}\) is the determinant. This formula has not been generalized to loop hafnians.

Tip

Implemented as thewalrus.hafnian() with approx=True. Note that one needs to pass the number of samples used to estimate the expectation value in the formula above; this is specified with the argument num_samples.

Low-rank algorithm

If a symmetric matrix \(\bm{A} \in \mathbb{C}^{n \times n}\) is of low rank it can be written as \(\bm{A} = \bm{G} \bm{G}^T\) where \(\bm{G} \in \mathbb{C}^{n \times r}\) is a rectangular matrix and \(r \leq n\) is the rank of \(\bm{A}\). One can then calculate the hafnian of the matrix \(\bm{A}\) in time and space complexity \({n+r-1 \choose r-1}\) by generalizing the result derived Barvinok [21] for permanents to hafnians as derived in Appendix C of Björklund et al [4]. The algorithm works by defining the following multivariate polynomial

\[q(x_0,\ldots,x_{r-1}) = \prod_{i=0}^{n-1} \sum_{j=0}^{r-1} G_{i,j} x_j.\]

Consider now the \(r\)-partitions of the integer \(n\), there are \({n+r-1 \choose r-1}\) of such partitions. Let \(e=(e_0,\ldots,e_{r-1})\) be an even \(r\)-partition (i.e. an \(r\)-partition where each element is even), call \(\mathcal{E}_{n,r}\) the set of all even \(r\)-partitions, and define \(\lambda_e\) to be the coefficient of \(x_0^{e_0}\ldots x_{r-1}^{e_{r-1}}\) in the polynomial \(q(x_0,\ldots,x_{r-1})\). Then one can write the hafnian of \(\bm{A} = \bm{G} \bm{G}^T\) as

\[\text{haf}(\bm{A}) = \sum_{e \in \mathcal{E}_{n,r}} \lambda_e \prod_{i=0}^{r-1} (e_i - 1)!! .\]

Tip

Implemented as thewalrus.low_rank_hafnian(). This function takes as argument the matrix \(\bm{G} \in \mathbb{C}^{n \times r}\) and returns the value of the hafnian of the matrix \(\bm{A} = \bm{G} \bm{G}^T\).

Sparse and banded algorithms

These algorithms take advantage of the Laplace expansion of the hafnian to calculate hafnians of matrices with many zeros. These matrices can be either banded or sparse. The Laplace expansion of the hafnian is given by

\[\text{haf}(\bm{B}) = \sum_{j \neq i} B_{i,j} \text{haf}(\bm{B}_{-i-j}),\]

where \(j\) is a fixed index and \(\bm{B}_{-i-j}\) is the matrix obtained from \(\bm{B}\) by removing rows and columns \(i\) and \(j\). The banded algorithm was introduced in Sec. V of Qi et al [22].

Tip

Implemented as thewalrus.hafnian_sparse() and thewalrus.hafnian_banded(). The loop hafnian calculation can be done by setting the option loop=True.

Multidimensional Hermite polynomials

Section author: Nicolás Quesada <nicolas@xanadu.ai>

In this section we study the multidimensional Hermite polynomials originally introduced by C. Hermite in 1865. See Mizrahi [23], Berkowitz et al. [24] and Kok and Braunstein [25] for more details.

In the next section, where we discuss quantum Gaussian states, we will explain how these polynomials relate to hafnians and loop hafnians. For the moment just let us introduce them and study their formal properties.

Generating function

Given two complex vectors \(\alpha,\beta \in \mathbb{C}^\ell\) and a symmetric matrix \(\bm{B} = \bm{B}^T \in \mathbb{C}^{\ell \times \ell}\),

\[\begin{split}F_B(\alpha,\beta) &= \exp\left( \alpha \bm{B} \beta^T - \tfrac{1}{2}\beta \bm{B} \beta^T\right) = \sum_{\bm{m} \geq \bm{0}} \prod_{i=1}^{\ell} \frac{\beta_i^{m_i}}{m_i!} H_{\bm{m}}^{(\bm{B})}(\alpha),\\ K_B(\alpha,\beta) &= \exp\left( \alpha \beta^T - \tfrac{1}{2}\beta \bm{B} \beta^T\right) = \sum_{\bm{m} \geq \bm{0}} \prod_{i=1}^{\ell} \frac{\beta_i^{m_i}}{m_i!} G_{\bm{m}}^{(\bm{B})}(\alpha),\end{split}\]

where the notation \(\bm{m} \geq \bm{0}\) is used to indicate that the sum goes over all vectors in \(\mathbb{N}^{\ell}_0\) (the set of vectors of nonnegative integers of size \(\ell\)). This generating function provides an implicit definition of the multidimensional Hermite polynomials. It is also straightforward to verify that \(H_{\bm{0}}^{(\bm{B})}(\alpha) = G_{\bm{0}}^{(\bm{B})}(\alpha) 1\). Finally, one can connect the standard Hermite polynomials \(H_{\bm{m}}^{(\bm{B})}(\alpha)\) to the modified Hermite polynomials \(G_{\bm{m}}^{(\bm{B})}(\alpha)\) via

\[H_{\bm{m}}^{(\bm{B})}(\alpha) = G_{\bm{m}}^{(\bm{B})}(\alpha \bm{B}).\]

In the one dimensional case, \(\ell=1\), one can compare the generating function above with the ones for the “probabilists’ Hermite polynomials” \(He_n(x)\) and “physicists’ Hermite polynomials” \(H_n(x)\) to find

\[\begin{split}He_n(x) = H_{n}^{([1])}(x), \\ H_n(x) = H_{n}^{([2])}(x).\end{split}\]

Tip

The standard multidimensional Hermite polynomials are implemented as thewalrus.hermite_multidimensional(). The modified Hermite polynomials can be obtained by passing the extra argument modified=True.

Recursion relation

Based on the generating functions introduced in the previous section one can derive the following recursion relations

\[\begin{split}H_{\bm{m}+\bm{e}_i}^{(\bm{B})}(\alpha) - \left(\sum_{j=1}^\ell B_{i,j} \alpha_j \right) H_{\bm{m}}^{(\bm{B})}(\alpha) + \sum_{j=1}^\ell B_{i,j} m_j H_{\bm{m}-\bm{e}_j}^{(\bm{B})}(\alpha) &= 0,\\ G_{\bm{m}+\bm{e}_i}^{(\bm{B})}(\alpha) - \alpha_i G_{\bm{m}}^{(\bm{B})}(\alpha) + \sum_{j=1}^\ell B_{i,j} m_j G_{\bm{m}-\bm{e}_j}^{(\bm{B})}(\alpha) &= 0,\end{split}\]

where \(\bm{e}_j\) is a vector with zeros in all its entries except in the \(i^{\text{th}}\) entry where it has a one.

From this recursion relation, or by Taylor expanding the generating function, one easily finds

\[\begin{split}H_{\bm{e}_i}^{(\bm{B})}(\alpha) &= \sum_{j=1}^\ell B_{ij} \alpha_j,\\ G_{\bm{e}_i}^{(\bm{B})}(\alpha) &= \alpha_i.\end{split}\]

Using this recursion relation one can calculate all the multidimensional Hermite polynomials up to a given cutoff.

The connection between the multidimensional Hermite polynomials and pure Gaussian states was reported by Wolf [26], and later by Kramer, Moshinsky and Seligman [27]. This same connection was also pointed out by Doktorov, Malkin and Man’ko in the context of vibrational modes of molecules [28]. Furthermore, this connection was later generalized to mixed Gaussian states by Dodonov, Man’ko and Man’ko [29]. These matrix elements have the form

\[C \times \frac{H_{\bm{m}}^{(\bm{B})}(\alpha)}{\sqrt{\bm{m}!}} = C \times \frac{G_{\bm{m}}^{(\bm{B})}(\alpha \bm{B})}{\sqrt{\bm{m}!}}.\]

To obtain the standard or modified Hermite polynomials renormalized by the square root of the factorial of its index \(\sqrt{\bm{m}!}\) one can pass the optional argument renorm=True.

Multidimensional Hermite polynomials and hafnians

By connecting the results in page 815 of Dodonov et al. [29] with the results in page 546 of Kan [6] one obtains the following relation between the hafnian and the multidimensional Hermite polynomials

\[H_{\bm{m}}^{(-\bm{B})}(\bm{0}) = G_{\bm{m}}^{(-\bm{B})}(\bm{0})= \text{haf}(\bm{B}_{\bm{m}}),\]

and moreover one can generalize it to

\[G_{\bm{m}}^{(-\bm{B})}(\alpha) = \text{lhaf}\left(\text{vid}(\bm{B}_{\bm{m}},\alpha_{\bm{m}})\right),\]

for loop hafnians. With these two identifications one can use the recursion relations of the multidimensional Hermite polynomials to calculate all the hafnians of the reductions of a given matrix up to a given cutoff.

With these observations and using the recursion relations for the Hermite polynomials and setting \(\bm{m}=\bm{1} - \bm{e}_i, \ \alpha = 0\) one easily derives the well known Laplace expansion for the hafnian (cf. Sec. 4.1 of [1])

\[\text{haf}(\bm{B}) = \sum_{j \neq i} B_{i,j} \text{haf}(\bm{B}_{-i-j}),\]

where \(j\) is a fixed index and \(\bm{B}_{-i-j}\) is the matrix obtained from \(\bm{B}\) by removing rows and columns \(i\) and \(j\).

Gaussian States in the Fock basis

Section author: Nicolás Quesada <nicolas@xanadu.ai>

In this section we show how the matrix elements of so-called Gaussian states in the Fock basis are related to the hafnians and loop hafnians introduced in the previous sections.

Note

The results presented in this page use the conventions introduced in the notation page and the reader is strongly encouraged to get familiar with them.

Wigner functions and Gaussian states

The quantum state \(\varrho\) of an \(\ell\)-mode system can be uniquely characterized by its Wigner function [30]

\[\label{Eq: Wigner} W(\vec \alpha; \varrho) = \int \frac{\text{d}\vec \xi}{\pi^{2\ell}} \text{Tr}[\varrho \hat D(\vec \xi)] \exp\left(\vec \alpha^T \bm{\Omega} \ \vec \xi\right),\]

where \(\vec \alpha = (\alpha_0,\ldots, \alpha_{\ell-1},\alpha_0^*,\ldots, \alpha_{\ell-1}^*)\) and similarly \(\vec \xi = (\xi_0,\ldots, \xi_{\ell-1},\xi_0^*,\ldots, \xi_{\ell-1}^*)\) are bivectors of complex amplitudes where the second half of the vector is the complex conjugate of the first half. The displacement operator is defined as \(\hat D(\xi):=\exp(\vec{\xi}^T \bm{\Omega} \hat \zeta)\), where \(\bm{\Omega}= \left[ \begin{smallmatrix} 0 & \mathbb{I} \\ -\mathbb{I} & 0 \end{smallmatrix} \right]\) is the symplectic form and \(\hat\zeta_j\) is an operator vector of the mode creation and annihilation operators. Recall, \(\ell\) being the number of modes, we have \(\hat\zeta_j=\hat a_j\) and \(\hat \zeta_{\ell+j}=\hat a_j^\dagger\) for \(j=1,\cdots,\ell\). These bosonic creation and annihilation operators satisfy the canonical commutation relations \([\hat a_i, \hat a_j]\) = 0 and \([\hat a_i, \hat a_j^\dagger] = \delta_{ij}\).

A quantum state is called Gaussian if its Wigner function is Gaussian [31]. Any multimode Gaussian state \(\varrho\) is completely parametrized by its first and second moments, namely the vector of means \(\vec{\beta}\) with components

\[\label{eq:WignerMeans} \vec \beta_j = \text{Tr}[\varrho\hat\zeta_j],\]

and the Wigner-covariance matrix \(\bm{\sigma}\) with entries

\[\label{eq:WignerCovMat} \bm{\sigma}_{jk} = \text{Tr}[\varrho \{\hat{\zeta}_j,\hat{\zeta}_k^\dagger \}]/2 - \vec \beta_j \vec \beta_k^*,\]

where \(\{x,y\} := xy +yx\) denotes the anti-commutator.

The variables \((\alpha_0,\ldots,\alpha_{\ell-1})\) are said to be complex normal distributed with mean \((\beta_0,\ldots,\beta_{\ell-1}) \in \mathbb{C}^{\ell}\) and covariance matrix \({\bm{\sigma}}\) [32] which furthermore needs to satisfy the uncertainty relation [33]

\[\label{eq:uncertainty} {\bm{\sigma}} + \bm{Z}/2 \geq 0,\]

where \(\bm{Z} = \left( \begin{smallmatrix} \mathbb{I} & 0\\ 0& -\mathbb{I} \end{smallmatrix} \right)\). The covariance matrix \(\bm{\sigma}\) is customarily parametrized in the following block form [32]

\[\begin{split}\bm{\sigma} = \left[\begin{array}{c|c} \bm{\Gamma} & \bm{C} \\ \hline \bm{C}^* & \bm{\Gamma}^* \end{array} \right],\end{split}\]

where \(\bm{\Gamma}\) is hermitian and positive definite, while \(\bm{C}\) is symmetric.

Gaussian states in the quadrature basis

Historically, Gaussian states are parametrized not in terms of the covariance matrix \(\bm{\sigma}\) of the complex amplitudes \(\alpha_j\), but rather in terms of its quadrature components, the canonical positions \(q_j\) and canonical momenta \(p_j\),

\[\alpha_j = \frac{1}{\sqrt{2 \hbar}} \left( q_j+ i p_j \right),\]

where \(\hbar\) is a positive constant. There are at least three conventions for the value of this constant, \(\hbar \in \{1/2,1,2 \}\).

It is convenient to write the following vector \(\bm{r} = (q_0,\ldots,q_{\ell-1},p_0,\ldots,p_{\ell-1})\) whose mean is related to the vector of means \(\vec \beta\) as

\[\begin{split}\vec \beta &= \frac{1}{\sqrt{\hbar}}\bm{R} \bm{r}, \\ \bm{\bar{r}} &= \sqrt{\hbar} \bm{R}^\dagger \vec \beta, \\ \bm{R} &= \frac{1}{\sqrt{2}}\begin{bmatrix} \mathbb{I} & i \mathbb{I} \\ \mathbb{I} & -i \mathbb{I} \end{bmatrix}.\end{split}\]

Similarly the complex normal covariance matrix \(\bm{\sigma}\) of the variables \((\alpha_0,\ldots,\alpha_{\ell-1})\) is related to the normal covariance matrix \(\bm{V}\) of the variables \(\bm{r} = (q_0,\ldots,q_{\ell-1},p_0,\ldots,p_{\ell-1})\) as

\[\begin{split}\bm{\sigma} &= \frac{1}{\hbar} \ \bm{R} \bm{V} \bm{R}^\dagger \\ \bm{V} &= {\hbar} \ \bm{R}^\dagger \bm{\sigma} \bm{R}.\end{split}\]

Tip

To convert between the complex covariance matrix \(\bm{\sigma}\) and the quadrature covariance matrix \(\bm{V}\) use the functions thewalrus.quantum.Qmat() and thewalrus.quantum.Covmat()

An important property of Gaussian states is that reduced (or marginal) states of a global Gaussian state are also Gaussian. This implies that the reduced covariance matrix of a subsystem of a Gaussian state together with a reduced vector of means fully characterize a reduced Gaussian state. The reduced covariance matrix for modes \(S = i_1,\ldots,i_n\) is obtained from the covariance matrix of the global state \(\bm{\sigma}\) or \(\bm{V}\) by keeping the columns and rows \(i_1,\ldots,i_n\) and \(i_1+\ell,\ldots,i_n+\ell\) of the original covariance matrix \(\bm{\sigma}\). Similarly one obtains the vector of means by keeping only entries \(i_1,\ldots,i_n\) and \(i_1+\ell,\ldots,i_n+\ell\) of the original vector of means \(\vec \beta\) or \(\bm{\bar{r}}\). Using the notation previously introduced, one can succinctly write the covariance matrix of modes \(S=i_1,\ldots,i_m\) as \(\bm{\sigma}_{(S)}\) or \(\bm{V}_{(S)}\) , and similarly the vector of means as \(\vec{\beta}_{(S)}\) or \(\bm{\bar{r}}_{(S)}\).

Tip

To obtain the reduced covariance matrix and vector of means for a certain subset of the modes use thewalrus.quantum.reduced_gaussian().

Note that for \(\bm{V}\) to be a valid quantum covariance matrix it needs to be symmetric and satisfy the uncertainty relation

\[\bm{V} + i \frac{\hbar}{2} \bm{\Omega} \geq 0.\]

Tip

To verify if a given quadrature covariance matrix is a valid quantum covariance matrix use the function thewalrus.quantum.is_valid_cov()

A Gaussian state is pure \(\varrho = \ket{\psi} \bra{\psi}\) if and only if \(\text{det}\left( \tfrac{2}{\hbar} \bm{V} \right) = 1\).

Tip

To verify if a given quadrature covariance matrix is a valid quantum covariance matrix and corresponds to a pure state use the function thewalrus.quantum.is_pure_cov()

Finally, there is a special subset of Gaussian states called classical whose covariance matrix satisfies

\[\bm{V} \geq \tfrac{\hbar}{2}\mathbb{I}.\]

This terminology is explained in the next section when sampling is discussed.

Tip

To verify if a given quadrature covariance matrix is a valid quantum covariance matrix and corresponds to a classical state use the function thewalrus.quantum.is_classical_cov()

Gaussian states in the Fock basis

In this section we use a generalization [34, 35] of the results of Hamilton et al. [36] by providing an explicit expression for Fock basis matrix elements \(\langle \bm{m} | \rho | \bm{n} \rangle\), \(\bm{n} = (n_0,\ldots, n_{\ell-1}), \bm{m} = (m_0,\ldots, m_{\ell-1})\), of an \(\ell\)-mode Gaussian state \(\rho\) with covariance matrix \(\bm{\sigma}\) and displacement vector \(\vec \beta\). Note that these matrix elements can also be calculated using multidimensional Hermite polynomials as shown by Dodonov et al. [29]. Depending on how many of these elements are required one can prefer to calculate loop hafnians or multidimensional Hermite polynomials. In particular if one only needs a few matrix elements it is more advantageous to use the formulas derived below. On the other hand if one requires all the matrix elements up to a certain Fock occupation cutoff it is more efficient to use the methods of Dodonov et al., which are also implemented in this library.

We first define the following useful quantities:

\[\begin{split}\bm{X} &= \begin{bmatrix} 0 & \mathbb{I} \\ \mathbb{I} & 0 \end{bmatrix} , \\ \bm{\Sigma} &= \bm{\sigma} +\tfrac{1}{2} \mathbb{I}_{2\ell},\\ T &=\frac{\exp\left(-\tfrac{1}{2} \vec \beta^\dagger \bm{\Sigma}^{-1} \vec \beta \right)}{ \sqrt{\text{det}(\bm{\Sigma}) \prod_{s=1}^\ell n_s! m_s!}},\\ \bm{p} &= (n_0,\ldots,n_{\ell-1},m_0,\ldots,m_{\ell-1}).\end{split}\]

We refer to \(\bm{\Sigma}\) as the Husimi covariance matrix.

As shown in detail in Appendix A of Ref. [35], the Fock matrix elements of a Gaussian state \(\rho\) are given by the expression

\[\label{Eq: lhaf} \langle \bm{m} | \rho | \bm{n} \rangle = T \times \text{lhaf}( \text{vid}(\bm{A}_{\bm{p}}, \gamma_{ \bm{p}}) ),\]

where \(\text{lhaf}\) is the loop hafnian and \(\text{vid}\) is the vector in diagonal notation introduced in the notation section.

Note that one can also obtain the probability of detecting a certain photon number pattern \(\bm{n} = (n_0,\ldots,n_{\ell-1})\) by calculating

\[p(\bm{n}|\varrho) = \langle \bm{n} | \varrho | \bm{n} \rangle.\]

Tip

To obtain the matrix element of gaussian state with quadrature covariance matrix \(\bm{V}\) and vector of means \(\bm{r}\) use the function thewalrus.quantum.density_matrix_element().

Tip

To obtain the Fock space density matrix of gaussian state with quadrature covariance matrix \(\bm{V}\) and vector of means \(\bm{r}\) use the function thewalrus.quantum.density_matrix().

In the case where the Gaussian state \(\varrho = |\psi \rangle \langle \psi|\) is pure then the matrix element

\[\langle \bm{m} | \varrho | \bm{n} \rangle = \langle \bm{m} | \psi \rangle \langle \psi| \bm{n} \rangle\]

factorizes into a product of two amplitudes. In Ref. [34] it was shown that the Fock amplitude of a gaussian state is also given by a loop hafnian. Then, for pure states the matrix \(\bm{\bar{A}} = \bm{\bar{B}} \oplus \bm{\bar{B}}^*\).

Tip

To obtain the overlap of a pure gaussian state with quadrature covariance matrix \(\bm{V}\) and vector of means \(\bm{r}\) and a given Fock state \(\langle \bm{n}|\) use the function thewalrus.quantum.pure_state_amplitude().

Tip

To obtain the Fock space state vector (ket) of a pure gaussian state with quadrature covariance matrix \(\bm{V}\) and vector of means \(\bm{r}\) use the function thewalrus.quantum.state_vector().

Gaussian states and threshold detection

In the last section we sketched how to obtain the probability that a certain photon-number outcome is obtained when a Gaussian state is measured with photon-number detectors. In this section we show how to obtain the analogous probability for the case of threshold detectors. These binary outcome detectors can only distinguish between the the vacuum state and occupied states, and thus for a single mode they are described by the POVM elements

\[\hat{\Pi}_0^{(i)} = \ket{0_i} \bra{0_i} \text{ and } \hat{\Pi}_1^{(i)} = 1_i - \hat{\Pi}_0^{(i)},\]

where \(\ket{0_i}\) is the vacuum state of mode \(i\) and \(1_i\) is the identity in the Hilbert space of mode \(i\).

For an \(\ell\) mode Gaussian state with zero mean, the outcome of threshold detection in all of its modes is described by a bitstring vector \(\bm{n} = (n_0,\ldots,n_{\ell-1})\) and the probability of the event is given by Born’s rule according to

\[\begin{split}p(\bm{n}|\varrho) &= \text{Tr} \left( \prod_{i=1}^{\ell} \Pi_{n_i}^{(i)} \varrho \right) = \frac{\text{tor} \left(\bm{O}_{\{ \bm{n}\}} \right)}{\sqrt{\text{det}(\bm{\sigma})}}, \\ \bm{O} &= \left(\mathbb{I}_{2\ell} - \bm{\Sigma}^{-1} \right)\end{split}\]

where \(\text{tor}\) is the Torontonian. For \(2 \ell \times 2 \ell\) matrix \(\bm{O}\) the Torontonian is defined as

\[\text{ltor}(\bm{O}) = \sum_{S \in P([\ell])} (-1)^{|S|} \frac{1}{\sqrt{\det\left(\mathbb{I} - \bm{O}_{(S)}\right)}}\]

The torontonian can be thought of as a generating function for hafnians (cf. the trace algorithm formula in algorithms section). The torontonian algorithm can be specified recursively to improve its performance [37].

The loop Torontonian is defined as

\[\text{ltor}(\bm{O}, \bm{\gamma}) = \sum_{S \in P([\ell])} (-1)^{|S|} \frac{exp \{\frac{1}{2} \gamma_{(S)} \det\left(\mathbb{I} - \bm{O}_{(S)}\right)^{-1} \gamma_{(S)}^* \}}{\sqrt{\det\left(\mathbb{I} - \bm{O}_{(S)}\right)}}\]

where \(\text{ltor}\) is the loop Torontonian, and \(\gamma=(\Sigma^{-1}\alpha)^*\) where \(\alpha\) is a \(2 \ell \times 2 \ell\) vector :cite: bulmer2022threshold.

Tip

The torontonian and loop torontonian are implemented as thewalrus.tor() and thewalrus.ltor() respectively.

Gaussian Boson Sampling

Section author: Nicolas Quesada <nicolas@xanadu.ai>

What is Gaussian Boson Sampling

Gaussian Boson Sampling was introduced in Ref. [36] as a problem that could potentially show how a non-universal quantum computer shows exponential speed ups over a classical computer. In the ideal scenario a certain number \(N\) of squeezed states are sent into \(M \times M\) interferometer and are then proved by photon-number resolving detectors. Hamilton et al. argue that under certain conjectures and assumptions it should be extremely inefficient for a classical computer to generate the samples that the quantum optical experiment just sketched generates by construction. Note that the setup described by Hamilton et al. is a special case of a Gaussian state. We consider Gaussian states, that except for having zero mean, are arbitrary.

In this section we describe a classical algorithm introduced in Ref. [38], that generates GBS samples in exponential time in the number of photons measured. This algorithm reduces the generation of a sample to the calculation of a chain of conditional probabilities, in which subsequent modes are interrogated as to how many photons should be detected conditioned on the detection record of previously interrogated modes. The exponential complexity of the algorithm stems from the fact that the conditional probabilities necessary to calculate a sample are proportional to hafnians of matrices of size \(2N\times 2N\), where \(N\) is the number of photons detected.

An exponential-time classical algorithm for GBS

As explained in the previous section the reduced or marginal density matrices of a Gaussian state are also Gaussian, and it is straightforward to compute their marginal covariance matrices given the covariance matrix of the global quantum state.

Let \(\bm{\Sigma}\) be the Husimi covariance matrix of the Gaussian state being measured with photon-number resolving detectors. Let \(S = (i_0,i_1,\ldots,i_{m-1})\) be a set of indices specifying a subset of the modes. In particular we write \(S=[k] = (0,1,2,\ldots, k-1)\) for the first \(k\) modes. We write \(\bm{\Sigma}_{(S)}\) to indicate the covariance matrix of the modes specified by the index set \(S\) and define \(\bm{A}^{S} := \bm{X} \left[\mathbb{I} - \left( \bm{\Sigma}_{(S)}\right)^{-1} \right]\).

As shown in the previous section, these matrices can be used to calculate photon-number probabilities as

\[p(\bm{N} = \bm{n}) = \frac{\text{haf}(\bm{A}^{S}_{(\bm{n})})}{ \bm{n}! \sqrt{\det(\bm{\Sigma}_{(S)})}}\]

where \(\bm{N}=\left\{N_{i_1},N_{i_2},\ldots,N_{i_m} \right\}\) is a random variable denoting the measurement outcome, and \(\bm{n} = \left\{n_{i_1},n_{i_2},\ldots,n_{i_m} \right\}\) is a set of integers that represent the actual outcomes of the photon number measurements and \(\bm{n}! = \prod_{j=0}^{m-1} n_j!\)

Now, we want to introduce an algorithm to generate samples of the random variable \(\{N_0,\ldots,N_{M-1}\}\) distributed according to the GBS probability distribution. To generate samples, we proceed as follows: First, we can always calculate the following probabilities

\[p(N_0=n_0) = \frac{\text{Haf}\left(\bm{A}^{[0]}_{(n_0)}\right)}{ n_0! \sqrt{\det(\bm{\Sigma}_{(0)})}},\]

for \(n_0=0,1,\ldots, n_{\max}\), where \(n_{\max}\) is a cut-off on the maximum number of photons we can hope to detect in this mode. Having constructed said probabilities, we can always generate a sample for the number of photons in the first mode. This will fix \(N_0 = n_0\). Now we want to sample from \(N_1|N_0=n_0\). To this end, we use the definition of conditional probability

\[p(N_1=n_1|N_0=n_0)= \frac{p(N_1=n_1,N_0=n_0)}{p(N_0=n_0)} =\frac{\text{haf}\left(\bm{A}^{[1]}_{(n_0,n_1)}\right)}{n_0! n_1! \sqrt{\det(\bm{\Sigma}_{([2])})}} \frac{1}{p(N_0=n_0)}.\]

We can, as before, use this formula to sample the number of photons in the second mode conditioned on observing \(n_0\) photons in the first mode. Note that the factor \(p(N_1=n_1)\) is already known from the previous step. By induction, the conditional probability of observing \(n_k\) photons in the \(k\)-th mode satisfies given observations of \(n_0,\ldots,n_{k-1}\) photons in the previous \(k\) modes is given by

\[\begin{split}p(N_k=n_k|N_{k-1}=n_{k-1},\ldots,N_0=n_0) &= \frac{p(N_k=n_k,N_{k-1}=n_{k-1},\ldots,N_0=n_0) }{p(N_{k-1}=n_{k-1},\ldots,N_0=n_0)} \\ &=\frac{\text{haf}\left(\bm{A}^{[k+1]}_{(n_0,n_2,\ldots,n_k)}\right)}{n_0! n_1! \ldots n_{k+1}! \sqrt{\det(\bm{\Sigma}_{([k])})}} \frac{1}{p(N_{k-1}=n_{k-1},\ldots,N_0=n_0)},\end{split}\]

where \(p(N_{k-1}=n_{k-1},\ldots,N_0=n_0)\) has already been calculated from previous steps. The algorithm then proceeds as follows: for mode \(k\), we use the previous equation to sample the number of photons in that mode conditioned on the number of photons in the previous \(k\) modes. Repeating this sequentially for all modes produces the desired sample.

Tip

To generate samples from a gaussian state specified by a quadrature covariance matrix use thewalrus.samples.generate_hafnian_sample().

Note that the above algorithm can also be generalized to states with finite means for which one only needs to provide the mean with the optional argument mean.

Threshold detection samples

Note the arguments presented in the previous section can also be generalized to threshold detection. In this case one simple need to replace \(\text{haf} \to \text{tor}\) and \(\bm{A}^{[k+1]}_{(n_0,n_2,\ldots,n_k)} \to \bm{O}^{[k+1]}_{(n_0,n_2,\ldots,n_k)}\) where \(\bm{O}^{S} = \left[\mathbb{I} - \left( \bm{\Sigma}_{(S)}\right)^{-1} \right]\).

Tip

To generate threshold samples from a gaussian state specified by a quadrature covariance matrix use thewalrus.samples.generate_torontonian_sample().

Sampling of classical states

In the previous section it was mentioned that states whose covariance matrix satisfies \(\bm{V} \geq \frac{\hbar}{2}\mathbb{I}\) are termed classical. These designation is due to the fact that for these states it is possible to obtain a polynomial (cubic) time algorithm to generate photon number or threshold samples [39].

Tip

To generate photon number or threshold samples from a classical gaussian state specified by a quadrature covariance matrix use thewalrus.samples.hafnian_sample_classical_state() or 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 a positive semidefinite matrix, for details on how this is done cf. Ref. [40].

Tip

To generate photon number samples from a thermal state parametrized by a positive semidefinite real matrix use the module thewalrus.csamples().

Notation

Section author: Nicolás Quesada <nicolas@xanadu.ai>

Matrices and vectors

In this section we introduce the notation that is used in the rest of the documentation.

We deal with so-called bivectors in which the second half of their components is the complex conjugate of the first half. Thus if \(\bm{u} = (u_1,\ldots u_\ell) \in \mathbb{C}^{\ell}\) then \(\vec{\alpha} = (\bm{u},\bm{u}^*) = (u_1,\ldots,u_\ell,u_1^*,\ldots,u_\ell^*)\) is a bivector. We use uppercase letters for (multi)sets such as \(S = \{1,1,2\}\).

Following standard Python and C convention, we index starting from zero, thus the first \(\ell\) natural numbers are \([\ell]:=\{0,1,\ldots,\ell-1\}\). We define the \(\text{diag}\) function as follows: When acting on a vector \(\bm{u}\) it returns a square diagonal matrix \(\bm{u}\). When acting on a square matrix it returns a vector with its diagonal entries.

Finally, we define the vector in diagonal (\(\text{vid}\)) operation, that takes a matrix \(\bm{A}\) of size \(n \times n\) and a vector \(\bm{u}\) of size \(n\) and returns the matrix

\[\text{vid}(\bm{A},\bm{u}) = \bm{A} - \text{diag}(\text{diag}( \bm{A})) + \text{diag}(\bm{u}),\]

which is simply the matrix \(\bm{A}\) with the vector \(\bm{u}\) placed along its diagonal.

The reduction operation

It is very useful to have compact notation to deal with matrices that are constructed by removing or repeating rows and column of a given primitive matrix. Imagine for example a given \(4 \times 4\) matrix

\[\begin{split}\bm{A} = \left( \begin{array}{cccc} A_{0,0} & A_{0,1} & A_{0,2} & A_{0,3} \\ A_{1,0} & A_{1,1} & A_{1,2} & A_{1,3} \\ A_{2,0} & A_{2,1} & A_{2,2} & A_{2,3} \\ A_{3,0} & A_{3,1} & A_{3,2} & A_{3,3} \\ \end{array} \right),\end{split}\]

and that you want to construct the following matrix

\[\begin{split}\bm{A}'= \left( \begin{array}{c|ccc||cc} A_{0,0} & A_{0,1} & A_{0,1} & A_{0,1} & A_{0,3} & A_{0,3} \\ \hline A_{1,0} & A_{1,1} & A_{1,1} & A_{1,1} & A_{1,3} & A_{1,3} \\ A_{1,0} & A_{1,1} & A_{1,1} & A_{1,1} & A_{1,3} & A_{1,3} \\ A_{1,0} & A_{1,1} & A_{1,1} & A_{1,1} & A_{1,3} & A_{1,3} \\ \hline \hline A_{3,0} & A_{3,1} & A_{3,1} & A_{3,1} & A_{3,3} & A_{3,3} \\ A_{3,0} & A_{3,1} & A_{3,1} & A_{3,1} & A_{3,3} & A_{3,3} \\ \end{array} \right),\end{split}\]

where the first row and column have been kept as they were, the second row and column have been repeated 3 times, the third row and column have been eliminated and the fourth and the last row and column have been repeated twice. To specify the number of repetitions (or elimination) of a given row-column we simply specify a vector of integers where each value tells us the number of repetitions, and we use the value 0 to indicate that a given row-column is removed. Thus defining \(\bm{m}=(1,3,0,2)\) we find its reduction by \(\bm{m}\) to be precisely the matrix in the last equation

\[\bm{A}' = \bm{A}_{\bm{n}}.\]

One can also define the reduction operation on vectors. For instance if \(\bm{u}=(u_0,u_1,u_2,u_3)\) and \(\bm{m}=(1,3,0,2)\) then \(\bm{u}_\bm{n} = (u_0,u_1,u_1,u_1,u_3,u_3)\).

Tip

Reduction on block matrices

When dealing with Gaussian states one typically encounters \(2\ell \times 2 \ell\) block matrices of the following form

\[\begin{split}\bm{A} = \left(\begin{array}{c|c} \bm{C} & \bm{D} \\ \hline \bm{E} & \bm{F} \\ \end{array} \right),\end{split}\]

where \(\bm{C},\bm{D},\bm{E},\bm{F}\) are of size \(\ell \times \ell\). Now imagine that one applies the reduction operation by a vector \(\bm{n} \in \mathbb{N}^{\ell}\) to each of the blocks. We introduce the following notation

\[\begin{split}\bm{A}_{(\bm{n})} = \left(\begin{array}{c|c} \bm{C}_{\bm{n}} & \bm{D}_{\bm{n}} \\ \hline \bm{E}_{\bm{n}} & \bm{F}_{\bm{n}} \\ \end{array} \right),\end{split}\]

where we have used the notation \((\bm{n})\) with the round brackets \(()\) to indicate that the reduction is applied to the blocks of the matrix \(\bm{A}\).

Similarly to block matrices, one can also define a reduction operator for bivectors. Thus if \(\vec \beta = (u_0,u_1,u_2,u_3,u_0^*,u_1^*,u_2^*,u_3^*)\) and \(\bm{m}=(1,3,0,2)\), then

\[\vec \beta_{(\bm{n} ) } = (u_0,u_1,u_1,u_1,u_3,u_3,u_0^*,u_1^*,u_1^*,u_1^*,u_3^*,u_3^*).\]

The reduction operation in terms of sets

A different way of specifying how many times a given row and column must me repeated is by giving a set in which we simply list the columns to be repeated. Thus for example the reduction index vector \(\bm{n} = (1,3,0,2)\) can alternatively be given as the multiset \(S=\{0,1,1,1,3,3 \}\) where the element 0 appears once to indicate the first row and column is repeated once, the index 1 appears three times to indicate that this row and column are repeated three times, etcetera.

Similarly for matrices of even size for which the following partition makes sense

\[\begin{split}\bm{A} = \left(\begin{array}{c|c} \bm{C} & \bm{D} \\ \hline \bm{E} & \bm{F} \\ \end{array} \right),\end{split}\]

where \(\bm{A}\) is of size \(2\ell \times 2\ell\) and \(\bm{C},\bm{D},\bm{E},\bm{F}\) are of size \(\ell \times \ell\) we define

\[\begin{split}\bm{A}_{(S)} = \left(\begin{array}{c|c} \bm{C}_S & \bm{D}_S \\ \hline \bm{E}_S & \bm{F}_S \\ \end{array} \right).\end{split}\]

This implies that if the index \(i\) appears \(m_i\) times in \(S\) then the columns \(i\) and \(i+\ell\) of \(\bm{A}\) will be repeated \(m_i\) times in \(\bm{A}_S\).

For instance if

\[\begin{split}\bm{A} = \left( \begin{array}{ccc|ccc} A_{0,0} & A_{0,1} & A_{0,2} & A_{0,3} & A_{0,4} & A_{0,5} \\ A_{1,0} & A_{1,1} & A_{1,2} & A_{1,3} & A_{1,4} & A_{1,5} \\ A_{2,0} & A_{2,1} & A_{2,2} & A_{2,3} & A_{2,4} & A_{2,5} \\ \hline A_{3,0} & A_{3,1} & A_{3,2} & A_{3,3} & A_{3,4} & A_{3,5} \\ A_{4,0} & A_{4,1} & A_{4,2} & A_{4,3} & A_{4,4} & A_{4,5} \\ A_{5,0} & A_{5,1} & A_{5,2} & A_{5,3} & A_{5,4} & A_{5,5} \\ \end{array} \right),\end{split}\]

and \(S=\{0,0,2,2,2\}\) one finds

\[\begin{split}\bm{A}_{(S)} = \left( \begin{array}{cc|ccc|cc|ccc} A_{0,0} & A_{0,0} & A_{0,2} & A_{0,2} & A_{0,2} & A_{0,3} & A_{0,3} & A_{0,5} & A_{0,5} & A_{0,5} \\ A_{0,0} & A_{0,0} & A_{0,2} & A_{0,2} & A_{0,2} & A_{0,3} & A_{0,3} & A_{0,5} & A_{0,5} & A_{0,5} \\ \hline A_{2,0} & A_{2,0} & A_{2,2} & A_{2,2} & A_{2,2} & A_{2,3} & A_{2,3} & A_{2,5} & A_{2,5} & A_{2,5} \\ A_{2,0} & A_{2,0} & A_{2,2} & A_{2,2} & A_{2,2} & A_{2,3} & A_{2,3} & A_{2,5} & A_{2,5} & A_{2,5} \\ A_{2,0} & A_{2,0} & A_{2,2} & A_{2,2} & A_{2,2} & A_{2,3} & A_{2,3} & A_{2,5} & A_{2,5} & A_{2,5} \\ \hline A_{3,0} & A_{3,0} & A_{3,2} & A_{3,2} & A_{3,2} & A_{3,3} & A_{3,3} & A_{3,5} & A_{3,5} & A_{3,5} \\ A_{3,0} & A_{3,0} & A_{3,2} & A_{3,2} & A_{3,2} & A_{3,3} & A_{3,3} & A_{3,5} & A_{3,5} & A_{3,5} \\ \hline A_{5,0} & A_{5,0} & A_{5,2} & A_{5,2} & A_{5,2} & A_{5,3} & A_{5,3} & A_{5,5} & A_{5,5} & A_{5,5} \\ A_{5,0} & A_{5,0} & A_{5,2} & A_{5,2} & A_{5,2} & A_{5,3} & A_{5,3} & A_{5,5} & A_{5,5} & A_{5,5} \\ A_{5,0} & A_{5,0} & A_{5,2} & A_{5,2} & A_{5,2} & A_{5,3} & A_{5,3} & A_{5,5} & A_{5,5} & A_{5,5} \\ \end{array} \right).\end{split}\]

The notation also extends in a straightforward fashion for bivectors. For example \(\vec \beta = (u_0,u_1,u_2,u_3,u_0^*,u_1^*,u_2^*,u_3^*)\) and \(S=\{1,1,2\}\) then \(\vec \beta_{(S)} = (u_1,u_1,u_2,u_1^*,u_1^*,u_2^*)\).

This notation becomes handy when describing certain algorithms for the calculation of the hafnian and torontonian introduced in the rest of the documentation.

Combining reduction and vector in diagonal

Here we show some basic examples of how the reduction and vector in diagonal operations work together

Consider the following matrix

\[\begin{split}\Sigma = \left( \begin{array}{ccc|ccc} A_{0,0} & A_{0,1} & A_{0,2} & B_{0,0} & B_{0,1} & B_{0,2} \\ A_{1,0} & A_{1,1} & A_{1,2} & B_{1,0} & B_{1,1} & B_{1,2} \\ A_{2,0} & A_{2,1} & A_{2,2} & B_{2,0} & B_{2,1} & B_{2,2} \\ \hline C_{0,0} & C_{0,1} & C_{0,2} & D_{0,0} & D_{0,1} & D_{0,2} \\ C_{1,0} & C_{1,1} & C_{1,2} & D_{1,0} & D_{1,1} & D_{1,2} \\ C_{2,0} & C_{2,1} & C_{2,2} & D_{2,0} & D_{2,1} & D_{2,2} \\ \end{array} \right),\end{split}\]

and bivector \(\vec{\beta} = (\beta_0,\beta_1,\beta_2,\beta_0^*,\beta_1^*,\beta_2^*)\) and we are given the index vector \(\bm{u} = (1,0,3)\). Then we can calculate the following

\[\begin{split}\Sigma_{(\bm{u})} &= \left( \begin{array}{cccc|cccc} A_{0,0} & A_{0,2} & A_{0,2} & A_{0,2} & B_{0,0} & B_{0,2} & B_{0,2} & B_{0,2} \\ A_{2,0} & A_{2,2} & A_{2,2} & A_{2,2} & B_{2,0} & B_{2,2} & B_{2,2} & B_{2,2} \\ A_{2,0} & A_{2,2} & A_{2,2} & A_{2,2} & B_{2,0} & B_{2,2} & B_{2,2} & B_{2,2} \\ A_{2,0} & A_{2,2} & A_{2,2} & A_{2,2} & B_{2,0} & B_{2,2} & B_{2,2} & B_{2,2} \\ \hline C_{0,0} & C_{0,2} & C_{0,2} & C_{0,2} & D_{0,0} & D_{0,2} & D_{0,2} & D_{0,2} \\ C_{2,0} & C_{2,2} & C_{2,2} & C_{2,2} & D_{2,0} & D_{2,2} & D_{2,2} & D_{2,2} \\ C_{2,0} & C_{2,2} & C_{2,2} & C_{2,2} & D_{2,0} & D_{2,2} & D_{2,2} & D_{2,2} \\ C_{2,0} & C_{2,2} & C_{2,2} & C_{2,2} & D_{2,0} & D_{2,2} & D_{2,2} & D_{2,2} \\ \end{array} \right),\\ \vec \beta_{(\bm{u})} &= (\beta_0,\beta_2,\beta_2,\beta_2,\beta_0^*,\beta_2^*,\beta_2^*,\beta_2^*),\end{split}\]

and finally write

\[\begin{split}\text{vid}(\Sigma_{(\bm{u})},\vec \beta_{(\bm{u})})= \left( \begin{array}{cccc|cccc} \beta_{0} & A_{0,2} & A_{0,2} & A_{0,2} & B_{0,0} & B_{0,2} & B_{0,2} & B_{0,2} \\ A_{2,0} & \beta_{2} & A_{2,2} & A_{2,2} & B_{2,0} & B_{2,2} & B_{2,2} & B_{2,2} \\ A_{2,0} & A_{2,2} & \beta_{2} & A_{2,2} & B_{2,0} & B_{2,2} & B_{2,2} & B_{2,2} \\ A_{2,0} & A_{2,2} & A_{2,2} & \beta_{2} & B_{2,0} & B_{2,2} & B_{2,2} & B_{2,2} \\ \hline C_{0,0} & C_{0,2} & C_{0,2} & C_{0,2} & \beta_{0}^* & D_{0,2} & D_{0,2} & D_{0,2} \\ C_{2,0} & C_{2,2} & C_{2,2} & C_{2,2} & D_{2,0} & \beta_{2}^* & D_{2,2} & D_{2,2} \\ C_{2,0} & C_{2,2} & C_{2,2} & C_{2,2} & D_{2,0} & D_{2,2} & \beta_{2}^* & D_{2,2} \\ C_{2,0} & C_{2,2} & C_{2,2} & C_{2,2} & D_{2,0} & D_{2,2} & D_{2,2} & \beta_{2}^* \\ \end{array} \right).\end{split}\]

Note that because there are repetitions, the diagonal elements of the matrix \(\bm{A}\) appear off diagonal in \(\bm{A}_{(\bm{n})}\) and also in \(\text{vid}(\bm{A}_{(\bm{n})},\vec{\beta}_{\bm{n}})\).

One can ignore the block structure of the matrix \(A\) and bivector \(\vec{\beta}\), and treat them as 6 dimensional objects and use an index vector of the same length. If we now define \(\bm{p} = (1,0,3,0,2,2)\) one finds

\[\begin{split}\Sigma_{\bm{p}} &= \left( \begin{array}{cccccccc} A_{0,0} & A_{0,2} & A_{0,2} & A_{0,2} & B_{0,1} & B_{0,1} & B_{0,2} & B_{0,2} \\ A_{2,0} & A_{2,2} & A_{2,2} & A_{2,2} & B_{2,1} & B_{2,1} & B_{2,2} & B_{2,2} \\ A_{2,0} & A_{2,2} & A_{2,2} & A_{2,2} & B_{2,1} & B_{2,1} & B_{2,2} & B_{2,2} \\ A_{2,0} & A_{2,2} & A_{2,2} & A_{2,2} & B_{2,1} & B_{2,1} & B_{2,2} & B_{2,2} \\ C_{1,0} & C_{1,2} & C_{1,2} & C_{1,2} & D_{1,1} & D_{1,1} & D_{1,2} & D_{1,2} \\ C_{1,0} & C_{1,2} & C_{1,2} & C_{1,2} & D_{1,1} & D_{1,1} & D_{1,2} & D_{1,2} \\ C_{2,0} & C_{2,2} & C_{2,2} & C_{2,2} & D_{2,1} & D_{2,1} & D_{2,2} & D_{2,2} \\ C_{2,0} & C_{2,2} & C_{2,2} & C_{2,2} & D_{2,1} & D_{2,1} & D_{2,2} & D_{2,2} \\ \end{array} \right),\\ \vec{\beta}_{\bm{p}}&=(\beta_0,\beta_2,\beta_2,\beta_2,\beta_1^*,\beta_1^*,\beta_2^*,\beta_2^*).\end{split}\]

Note that we wrote \(\Sigma_{\bm{p}}\) and not \(\Sigma_{(\bm{p})}\) to indicate that we ignore the block structure of the matrix \(\Sigma\).

References

[1]

Alexander Barvinok. Combinatorics and complexity of partition functions. Volume 274. Springer, 2016.

[2]

Eduardo R Caianiello. On quantum field theory—i: explicit solution of dyson’s equation in electrodynamics without use of feynman graphs. Il Nuovo Cimento (1943-1954), 10(12):1634–1652, 1953.

[3]

Settimo Termini. Imagination and Rigor: Essays on Eduardo R. Caianiello's Scientific Heritage. Springer, 2006.

[4]

Andreas Björklund, Brajesh Gupt, and Nicolás Quesada. A faster hafnian formula for complex matrices and its benchmarking on a supercomputer. Journal of Experimental Algorithmics (JEA), 24(1):11, 2019.

[5]

Andreas Björklund and Thore Husfeldt. Exact algorithms for exact satisfiability and number of perfect matchings. Algorithmica, 52(2):226–249, 2008.

[6]

Raymond Kan. From moments of sum to moments of product. Journal of Multivariate Analysis, 99(3):542–554, 2008.

[7]

Mikko Koivisto. Partitioning into sets of bounded cardinality. In International Workshop on Parameterized and Exact Computation, 258–263. Springer, 2009.

[8]

Jesper Nederlof. Fast polynomial-space algorithms using möbius inversion: improving on steiner tree and related problems. In International Colloquium on Automata, Languages, and Programming, 713–725. Springer, 2009.

[9]

Andreas Björklund. Counting perfect matchings as fast as ryser. In Proceedings of the twenty-third annual ACM-SIAM symposium on Discrete Algorithms, 914–921. SIAM, 2012.

[10]

Marek Cygan and Marcin Pilipczuk. Faster exponential-time algorithms in graphs of bounded average degree. Information and Computation, 243:75–85, 2015.

[11]

Herbert John Ryser. Combinatorial mathematics. Number 14. Mathematical Association of America; distributed by Wiley [New York], 1963.

[12]

Grzegorz Rempala and Jacek Wesolowski. Symmetric functionals on random matrices and random matchings problems. Volume 147. Springer Science & Business Media, 2007.

[13]

Donald E. Knuth. The Art of Computer Programming, Vol. 2: Seminumerical Algorithms. Addison-Wesley, 3 edition, 1998.

[14]

Leslie G Valiant. The complexity of computing the permanent. Theoretical computer science, 8(2):189–201, 1979.

[15]

Alexander Barvinok. Approximating permanents and hafnians. arXiv preprint arXiv:1601.07518, 2016.

[16]

Piotr Sankowski. Alternative algorithms for counting all matchings in graphs. In Helmut Alt and Michel Habib, editors, STACS 2003, 427–438. Berlin, Heidelberg, 2003. Springer Berlin Heidelberg.

[17]

Alexander Barvinok. Polynomial time algorithms to approximate permanents and mixed discriminants within a simply exponential factor. Random Structures & Algorithms, 14(1):29–61, 1999.

[18]

Mark Rudelson, Alex Samorodnitsky, Ofer Zeitouni, and others. Hafnians, perfect matchings and gaussian matrices. The Annals of Probability, 44(4):2858–2888, 2016.

[19]

Nicolás Quesada, Juan Miguel Arrazola, and Nathan Killoran. Gaussian boson sampling using threshold detectors. Physical Review A, 98(6):062322, 2018.

[20]

Rizwana Rehman and Ilse CF Ipsen. La budde's method for computing characteristic polynomials. arXiv preprint arXiv:1104.3769, 2011.

[21]

Alexander I Barvinok. Two algorithmic results for the traveling salesman problem. Mathematics of Operations Research, 21(1):65–84, 1996.

[22]

Haoyu Qi, Diego Cifuentes, Kamil Brádler, Robert Israel, Timjan Kalajdzievski, and Nicolás Quesada. Efficient sampling from shallow gaussian quantum-optical circuits with local interactions. arXiv preprint arXiv:2009.11824, 2020.

[23]

Maurice M Mizrahi. Generalized hermite polynomials. Journal of Computational and Applied Mathematics, 1(3):137–140, 1975.

[24]

S Berkowitz and FJ Garner. The calculation of multidimensional Hermite polynomials and Gram-Charlier coefficients. Math. Comput., 24(111):537–545, 1970.

[25]

Pieter Kok and Samuel L Braunstein. Multi-dimensional hermite polynomials in quantum optics. Journal of Physics A: Mathematical and General, 34(31):6185, 2001.

[26]

Kurt Bernardo Wolf. Canonical transforms. i. complex linear transforms. Journal of Mathematical Physics, 15(8):1295–1301, 1974.

[27]

Kramer P, M Moshinsky, and T.H Seligman. Complex extensions of canonical transformations and quantum mechanics. In E.M. Loebl, editor, Group Theory and Its Applications, pages 250–300. Academic, New York, 1975.

[28]

EV Doktorov, IA Malkin, and VI Man'ko. Dynamical symmetry of vibronic transitions in polyatomic molecules and the franck-condon principle. Journal of Molecular Spectroscopy, 64(2):302–326, 1977.

[29]

VV Dodonov, OV Man’ko, and VI Man’ko. Multidimensional hermite polynomials and photon distribution for polymode mixed light. Physical Review A, 50(1):813, 1994.

[30]

Alessio Serafini. Quantum Continuous Variables: A Primer of Theoretical Methods. CRC Press, 2017.

[31]

Christian Weedbrook, Stefano Pirandola, Raúl García-Patrón, Nicolas J Cerf, Timothy C Ralph, Jeffrey H Shapiro, and Seth Lloyd. Gaussian quantum information. Reviews of Modern Physics, 84(2):621, 2012.

[32]

Bernard Picinbono. Second-order complex random vectors and normal distributions. IEEE Transactions on Signal Processing, 44(10):2637–2640, 1996.

[33]

R Simon, N Mukunda, and Biswadeb Dutta. Quantum-noise matrix for multimode systems: $U(n)$ invariance, squeezing, and normal forms. Physical Review A, 49(3):1567, 1994.

[34]

Nicolás Quesada. Franck-condon factors by counting perfect matchings of graphs with loops. The Journal of chemical physics, 150(16):164113, 2019.

[35]

N. Quesada, L. G. Helt, J. Izaac, J. M. Arrazola, R. Shahrokhshahi, C. R. Myers, and K. K. Sabapathy. Simulating realistic non-gaussian state preparation. Phys. Rev. A, 100:022341, 2019.

[36]

Craig S Hamilton, Regina Kruse, Linda Sansoni, Sonja Barkhofen, Christine Silberhorn, and Igor Jex. Gaussian boson sampling. Physical review letters, 119(17):170501, 2017.

[37]

Ágoston Kaposi, Zoltán Kolarovszki, Tamás Kozsik, Zoltán Zimborás, and Péter Rakyta. Polynomial speedup in torontonian calculation by a scalable recursive algorithm. arXiv preprint arXiv:2109.04528, 2021.

[38]

Nicolás Quesada and Juan Miguel Arrazola. Exact simulation of gaussian boson sampling in polynomial space and exponential time. Physical Review Research, 2(2):023005, 2020.

[39]

Saleh Rahimi-Keshari, Austin P Lund, and Timothy C Ralph. What can quantum optics say about computational complexity theory? Physical review letters, 114(6):060501, 2015.

[40]

Soran Jahangiri, Juan Miguel Arrazola, Nicolás Quesada, and Nathan Killoran. Point processes with gaussian boson sampling. Physical Review E, 101(2):022134, 2020.

[41]

Francesco Mezzadri. How to generate random matrices from the classical compact groups. arXiv e-prints, pages math–ph/0609050, September 2006. arXiv:math-ph/0609050.

[42]

Martin Houde, Will McCutcheon, and Nicolás Quesada. Matrix decompositions in quantum optics: takagi/autonne, bloch-messiah/euler, iwasawa, and williamson. arXiv preprint arXiv:2403.04596, 2024.

Overview

Python API

  • The thewalrus Python API provides access to various highly optimized hafnian, permanent, and torontonian algorithms

  • The thewalrus.quantum submodule provides access to various utility functions that act on Gaussian quantum states

  • The thewalrus.samples submodule provides access to algorithms to sample from the hafnian or the torontonian of Gaussian quantum states

  • The thewalrus.symplectic submodule provides access to a convenient set of symplectic matrices and utility functions to manipulate them

  • The thewalrus.decompositions submodule provides access to common shared matrix decompositions used to perform gate decompositions

  • The thewalrus.charpoly submodule provides access to La Budde’s algorithm to calculate the characteristic polynomials of matrices

  • The thewalrus.random submodule provides access to random unitary, symplectic and covariance matrices

  • The thewalrus.fock_gradients submodule provides access to the Fock representation of certain continuous-variable gates and their gradients

  • The thewalrus.reference submodule provides access to pure-Python reference implementations of the hafnian, loop hafnian, and torontonian

Octave

In addition, two Octave functions are provided: octave/hafnian.m and octave/loophafnian.m.

The Walrus

This is the top level module of the The Walrus, containing functions for computing the hafnian, loop hafnian, and torontonian of matrices.

Algorithm terminology

Eigenvalue hafnian algorithm

The algorithm described in A faster hafnian formula for complex matrices and its benchmarking on a supercomputer, [4]. This algorithm scales like \(\mathcal{O}(n^3 2^{n/2})\), and supports calculation of the loop hafnian.

Recursive hafnian algorithm

The algorithm described in Counting perfect matchings as fast as Ryser [9]. This algorithm scales like \(\mathcal{O}(n^4 2^{n/2})\). This algorithm does not currently support the loop hafnian.

Repeated hafnian algorithm

The algorithm described in From moments of sum to moments of product, [6]. This method is more efficient for matrices with repeated rows and columns, and supports calculation of the loop hafnian.

Approximate hafnian algorithm

The algorithm described in Polynomial time algorithms to approximate permanents and mixed discriminants within a simply exponential factor, [17]. This algorithm allows us to efficiently approximate the hafnian of matrices with non-negative elements. This is done by sampling determinants; the larger the number of samples taken, the higher the accuracy.

Batched hafnian algorithm

An algorithm that allows the calculation of hafnians of all reductions of a given matrix up to the cutoff (resolution) provided. Internally, this algorithm makes use of the multidimensional Hermite polynomials as per The calculation of multidimensional Hermite polynomials and Gram-Charlier coefficients [24].

Low-rank hafnian algorithm

An algorithm that allows to calculate the hafnian of an \(r\)-rank matrix \(\bm{A}\) of size \(n \times n\) by factorizing it as \(\bm{A} = \bm{G} \bm{G}^T\) where \(\bm{G}\) is of size \(n \times r\). The algorithm is described in Appendix C of A faster hafnian formula for complex matrices and its benchmarking on a supercomputer, [4].

Banded hafnian algorithm

An algorithm that calculates the hafnian of a matrix \(\bm{A}\) of size \(n \times n\) with bandwidth \(w\) by calculating and storing certain subhafnians dictated by the bandwidth. The algorithm is described in Section V of Efficient sampling from shallow Gaussian quantum-optical circuits with local interactions, [22].

Sparse hafnian algorithm

An algorithm that calculates the hafnian of a sparse matrix by taking advantage of the Laplace expansion and memoization, to store only the relevant paths that contribute non-zero values to the final calculation.

Functions

hafnian(A[, loop, rtol, atol, approx, ...])

Returns the hafnian of a matrix.

hafnian_repeated(A, rpt[, mu, loop, rtol, ...])

Returns the hafnian of matrix with repeated rows/columns.

hafnian_batched(A, cutoff[, mu, rtol, atol, ...])

Calculates the hafnian of reduction(A, k) for all possible values of vector k below the specified cutoff.

tor(A[, recursive])

Returns the Torontonian of a matrix.

ltor(A, gamma[, recursive])

Returns the loop Torontonian of an NxN matrix and an N-length vector.

perm(A[, method])

Returns the permanent of a matrix using various methods.

permanent_repeated(A, rpt)

Calculates the permanent of matrix \(A\), where the ith row/column of \(A\) is repeated \(rpt_i\) times.

brs(A, E)

Calculates the Bristolian, a matrix function introduced for calculating the threshold detector statistics on measurements of Fock states interfering in linear optical interferometers.

ubrs(A)

Calculates the Unitary Bristolian, a matrix function introduced for calculating the threshold detector statistics on measurements of Fock states interfering in lossless linear optical interferometers.

hermite_multidimensional(R, cutoff[, y, C, ...])

Returns photon number statistics of a Gaussian state for a given covariance matrix as described in Multidimensional Hermite polynomials and photon distribution for polymode mixed light arxiv:9308033.

hafnian_banded(A[, loop, rtol, atol])

Returns the loop hafnian of a banded matrix.

reduction(A, rpt)

Calculates the reduction of an array by a vector of indices.

version()

Get version number of The Walrus

low_rank_hafnian(G)

Returns the hafnian of the low rank matrix \(\bm{A} = \bm{G} \bm{G}^T\) where \(\bm{G}\) is rectangular of size \(n \times r\) with \(r \leq n\).

Code details

brs(A, E)[source]

Calculates the Bristolian, a matrix function introduced for calculating the threshold detector statistics on measurements of Fock states interfering in linear optical interferometers.

See the paper ‘Threshold detector statistics of Bosonic states’ for more detail (to be published soon)

Parameters:

  • A (array) – matrix of size [m, n]

  • E (array) – matrix of size [r, n]

Returns:

the Bristol of matrices A and E

Return type:

int or float or complex

grad_hermite_multidimensional(G, R, y, C=1, renorm=True, dtype=None)[source]

Calculates the gradients of the renormalized multidimensional Hermite polynomials \(C*H_k^{(R)}(y)\) with respect to its parameters \(C\), \(y\) and \(R\).

Parameters:

  • G (array) – the multidimensional Hermite polynomials

  • R (array[complex]) – square matrix parametrizing the Hermite polynomial

  • y (vector[complex]) – vector argument of the Hermite polynomial

  • C (complex) – first value of the Hermite polynomials

  • renorm (bool) – If True, uses the normalized multidimensional Hermite polynomials such that \(H_k^{(R)}(y)/\prod_i k_i!\)

  • dtype (data type) – Specifies the data type used for the calculation

Returns:

the gradients of the multidimensional Hermite polynomials with respect to C, R and y

Return type:

array[data type], array[data type], array[data type]

hafnian(A, loop=False, rtol=1e-05, atol=1e-08, approx=False, num_samples=1000, method='glynn')[source]

Returns the hafnian of a matrix. Code contributed by Jake F.F. Bulmer based on arXiv:2108.01622.

Parameters:

  • A (array) – a square, symmetric array of even dimensions

  • loop (bool) – If True, the loop hafnian is returned. Default is False.

  • method (string) – Set this to "glynn" to use the glynn formula, or "inclexcl" to use the inclusion exclusion principle, or "recursive" to use a recursive algorithm.

  • rtol (float) – the relative tolerance parameter used in np.allclose

  • atol (float) – the absolute tolerance parameter used in np.allclose

  • approx (bool) – If True, an approximation algorithm is used to estimate the hafnian. Note that the approximation algorithm can only be applied to matrices A that only have non-negative entries.

  • num_samples (int) – If approx=True, the approximation algorithm performs num_samples iterations for estimation of the hafnian of the non-negative matrix A

Returns:

the hafnian of matrix A

Return type:

int or float or complex

hafnian_banded(A, loop=False, rtol=1e-05, atol=1e-08)[source]

Returns the loop hafnian of a banded matrix. For the derivation see Section V of ‘Efficient sampling from shallow Gaussian quantum-optical circuits with local interactions’, Qi et al..

Parameters:

A (array) – a square, symmetric array of even dimensions

Returns:

the loop hafnian of matrix A

Return type:

int or float or complex

hafnian_batched(A, cutoff, mu=None, rtol=1e-05, atol=1e-08, renorm=False, make_tensor=True)[source]

Calculates the hafnian of reduction(A, k) for all possible values of vector k below the specified cutoff.

Here,

  • \(A\) is am \(n\times n\) square matrix

  • \(k\) is a vector of (non-negative) integers with the same dimensions as \(A\), i.e., \(k = (k_0,k_1,\ldots,k_{n-1})\), and where \(0 \leq k_j < \texttt{cutoff}\).

The function hafnian_repeated() can be used to calculate the reduced hafnian for a specific value of \(k\); see the documentation for more information.

Note

If mu is not None, this function instead returns hafnian(np.fill_diagonal(reduction(A, k), reduction(mu, k)), loop=True). This calculation can only be performed if the matrix \(A\) is invertible.

Parameters:

  • A (array) – a square, symmetric \(N\times N\) array.

  • cutoff (int) – maximum size of the subindices in the Hermite polynomial

  • mu (array) – a vector of length \(N\) representing the vector of means/displacement

  • renorm (bool) – If True, the returned hafnians are normalized, that is, \(haf(reduction(A, k))/\ \sqrt{prod_i k_i!}\) (or \(lhaf(fill\_diagonal(reduction(A, k),reduction(mu, k)))\) if mu is not None)

  • make_tensor – If False, returns a flattened one dimensional array instead of a tensor with the values of the hafnians.

  • rtol (float) – the relative tolerance parameter used in np.allclose.

  • atol (float) – the absolute tolerance parameter used in np.allclose.

Returns:

the values of the hafnians for each value of \(k\) up to the cutoff

Return type:

(array)

hafnian_repeated(A, rpt, mu=None, loop=False, rtol=1e-05, atol=1e-08, glynn=True)[source]

Returns the hafnian of matrix with repeated rows/columns.

Code contributed by Jake F.F. Bulmer based on arXiv:2108.01622.

The reduction() function may be used to show the resulting matrix with repeated rows and columns as per rpt.

As a result, the following are identical:

However, using hafnian_repeated in the case where there are a large number of repeated rows and columns (\(\sum_{i}rpt_i \gg N\)) can be significantly faster.

Note

If \(rpt=(1, 1, \dots, 1)\), then

>>> hafnian_repeated(A, rpt) == hafnian(A)
Parameters:

  • A (array) – a square, symmetric \(N\times N\) array

  • rpt (Sequence) – a length-\(N\) positive integer sequence, corresponding to the number of times each row/column of matrix \(A\) is repeated

  • mu (array) – A vector of length \(N\) representing the vector of means/displacement. If not provided, mu is set to the diagonal of matrix A. Note that this only affects the loop hafnian.

  • loop (bool) – If True, the loop hafnian is returned. Default is False.

  • rtol (float) – the relative tolerance parameter used in np.allclose

  • atol (float) – the absolute tolerance parameter used in np.allclose

  • glynn (bool) – whether to use finite difference sieve

Returns:

the hafnian of matrix A

Return type:

int or float or complex

hafnian_sparse(A, D=None, loop=False)[source]

Returns the hafnian of a sparse symmetric matrix. This pure python implementation is very slow on full matrices, but faster the sparser a matrix is. As a rule of thumb, the crossover in runtime with respect to hafnian() happens around 50% sparsity.

Parameters:

  • A (array) – the symmetric matrix of which we want to compute the hafnian

  • D (set) – Set of indices that identify a submatrix. If None (default) it computes the hafnian of the whole matrix.

  • loop (bool) – If True, the loop hafnian is returned. Default is False.

Returns:

hafnian of A or of the submatrix of A defined by the set of indices D

Return type:

float

hermite_multidimensional(R, cutoff, y=None, C=1, renorm=False, make_tensor=True, modified=False, rtol=1e-05, atol=1e-08)[source]

Returns photon number statistics of a Gaussian state for a given covariance matrix as described in Multidimensional Hermite polynomials and photon distribution for polymode mixed light arxiv:9308033.

Here \(R\) is an \(n \times n\) square matrix, and \(y\) is an \(n\) dimensional vector. The polynomials \(H_k^{(R)}(y)\) are parametrized by the multi-index \(k=(k_0,k_1,\ldots,k_{n-1})\), and are calculated for all values \(0 \leq k_j < \text{cutoff}\), thus a tensor of dimensions \(\text{cutoff}^n\) is returned.

This tensor can either be flattened into a vector or returned as an actual tensor with \(n\) indices.

This implementation is based on the MATLAB code available at github clementsw/gaussian-optics.

Note

Note that if \(R = (1)\) then \(H_k^{(R)}(y)\) are precisely the well known probabilists’ Hermite polynomials \(He_k(y)\), and if \(R = (2)\) then \(H_k^{(R)}(y)\) are precisely the well known physicists’ Hermite polynomials \(H_k(y)\).

Parameters:

  • R (array) – square matrix parametrizing the Hermite polynomial family

  • cutoff (int) – maximum size of the subindices in the Hermite polynomial

  • y (array) – vector argument of the Hermite polynomial

  • C (complex) – first value of the Hermite polynomials, the default value is 1

  • renorm (bool) – If True, normalizes the returned multidimensional Hermite polynomials such that \(H_k^{(R)}(y)/\prod_i k_i!\)

  • make_tensor (bool) – If False, returns a flattened one dimensional array containing the values of the polynomial

  • modified (bool) – whether to return the modified multidimensional Hermite polynomials or the standard ones

  • rtol (float) – the relative tolerance parameter used in np.allclose

  • atol (float) – the absolute tolerance parameter used in np.allclose

Returns:

the multidimensional Hermite polynomials

Return type:

(array)

lmtl(A, zeta)[source]

Returns the montrealer of an NxN matrix and an N-length vector.

Parameters:

  • A (array) – an NxN array of even dimensions

  • zeta (array) – an N-length vector of even dimensions

Returns:

the loop montrealer of matrix A, vector zeta

Return type:

np.float64 or np.complex128

loop_hafnian(A, D=None, reps=None, glynn=True)[source]

Calculate loop hafnian with (optional) repeated rows and columns. Code contributed by Jake F.F. Bulmer based on arXiv:2108.01622.

Parameters:

  • A (array) – N x N matrix.

  • D (array) – Diagonal entries of matrix (optional). If not provided, D is the diagonal of A. If repetitions are provided, D should be provided explicitly.

  • reps (list) – Length-N list of repetitions of each row/col (optional), if not provided, each row/column assumed to be repeated once.

  • glynn (bool) – If True, use Glynn-style finite difference sieve formula, if False, use Ryser style inclusion/exclusion principle.

Returns

complex: result of loop hafnian calculation

ltor(A, gamma, recursive=True)[source]

Returns the loop Torontonian of an NxN matrix and an N-length vector.

Parameters:

  • A (array) – an NxN array of even dimensions.

  • gamma (array) – an N-length vector of even dimensions

  • recursive – use the faster recursive implementation

Returns:

the loop torontonian of matrix A, vector gamma

Return type:

np.float64 or np.complex128

matched_reps(reps)[source]

Takes the repeated rows and find a way to pair them up to create a perfect matching with many repeated edges.

Parameters:

reps (list) – list of repeated rows/cols

Returns:

tuple with vertex pairs (length 2N for N edges; index i is matched with i + N), length N array for how many times each edge is repeated, and index of odd mode (None if even number of vertices)

Return type:

tuple[array, array, int]

mtl(A)[source]

Returns the montrealer of an NxN matrix.

Parameters:

A (array) – an NxN array of even dimensions.

Returns:

the montrealer of matrix A

Return type:

np.float64 or np.complex128

perm(A, method='bbfg')[source]

Returns the permanent of a matrix using various methods.

Parameters:

  • A (array[float or complex]) – a square array.

  • method (string) – Set this to "ryser" to use the Ryser formula, or "bbfg" to use the BBFG formula.

Returns:

the permanent of matrix A

Return type:

float or complex

permanent_repeated(A, rpt)[source]

Calculates the permanent of matrix \(A\), where the ith row/column of \(A\) is repeated \(rpt_i\) times.

This function constructs the matrix

\[\begin{split}B = \begin{bmatrix} 0 & A\\ A^T & 0 \end{bmatrix},\end{split}\]

and then calculates \(perm(A)=haf(B)\), by calling

>>> hafnian_repeated(B, rpt*2, loop=False)
Parameters:

  • A (array) – matrix of size [N, N]

  • rpt (Sequence) – sequence of N positive integers indicating the corresponding rows/columns of A to be repeated.

Returns:

the permanent of matrix A

Return type:

int or float or complex

reduction(A, rpt)[source]

Calculates the reduction of an array by a vector of indices. This is equivalent to repeating the ith row/column of \(A\), \(rpt_i\) times.

Parameters:

  • A (array) – matrix of size [N, N]

  • rpt (Sequence) – sequence of N positive integers indicating the corresponding rows/columns of A to be repeated.

Returns:

the reduction of A by the index vector rpt

Return type:

array

tor(A, recursive=True)[source]

Returns the Torontonian of a matrix.

Parameters:

  • A (array) – a square array of even dimensions.

  • recursive – use the faster recursive implementation.

Returns:

the torontonian of matrix A.

Return type:

np.float64 or np.complex128

ubrs(A)[source]

Calculates the Unitary Bristolian, a matrix function introduced for calculating the threshold detector statistics on measurements of Fock states interfering in lossless linear optical interferometers.

See the paper ‘Threshold detector statistics of Bosonic states’ for more detail (to be published soon)

Parameters:

A (array) – matrix of size [m, n]

Returns:

the Unitary Bristol of matrix A

Return type:

int or float or complex

version()[source]

Get version number of The Walrus

Returns:

The package version number

Return type:

str

Quantum algorithms

Module name: thewalrus.quantum

This submodule provides access to various utility functions that act on Gaussian quantum states.

For more details on how the hafnian relates to various properties of Gaussian quantum states, see:

Fock states and tensors

pure_state_amplitude(mu, cov, i[, ...])

Returns the \(\langle i | \psi\rangle\) element of the state ket of a Gaussian state defined by covariance matrix cov.

state_vector(mu, cov[, post_select, ...])

Returns the state vector of a (PNR post-selected) Gaussian state.

density_matrix_element(mu, cov, i, j[, ...])

Returns the \(\langle i | \rho | j \rangle\) element of the density matrix of a Gaussian state defined by covariance matrix cov.

density_matrix(mu, cov[, post_select, ...])

Returns the density matrix of a (PNR post-selected) Gaussian state.

fock_tensor(S, alpha, cutoff[, choi_r, ...])

Calculates the Fock representation of a Gaussian unitary parametrized by the symplectic matrix S and the displacements alpha up to cutoff in Fock space.

probabilities(mu, cov, cutoff[, parallel, ...])

Generate the Fock space probabilities of a Gaussian state up to a Fock space cutoff.

loss_mat(eta, cutoff)

Constructs a binomial loss matrix with transmission eta up to n photons.

update_probabilities_with_loss(etas, probs)

Given a list of transmissivities a tensor of probabilitites, calculate an updated tensor of probabilities after loss is applied.

update_probabilities_with_noise(probs_noise, ...)

Given a list of noise probability distributions for each of the modes and a tensor of probabilitites, calculate an updated tensor of probabilities after noise is applied.

find_classical_subsystem(cov[, hbar, atol])

Find the largest integer k so that subsystem in modes [0,1,...,k-1] is a classical state.

tvd_cutoff_bounds(mu, cov, cutoff[, hbar, ...])

Gives bounds of the total variation distance between the exact Gaussian Boson Sampling distribution extending to infinity in Fock space and the ones truncated by any value between 0 and the user provided cutoff.

n_body_marginals(mean, cov, cutoff, n[, hbar])

Calculates the first n-body marginals of a Gaussian state.

Adjacency matrices

adj_scaling(A, n_mean)

Returns the scaling parameter by which the adjacency matrix A should be rescaled so that the Gaussian state that endodes it has a total mean photon number n_mean.

adj_scaling_torontonian(A, c_mean)

Returns the scaling parameter by which the adjacency matrix A should be rescaled so that the Gaussian state that encodes it has give a mean number of clicks equal to c_mean when measured with threshold detectors.

adj_to_qmat(A, n_mean)

Returns the Qmat xp-covariance matrix associated to a graph with adjacency matrix \(A\) and with mean photon number \(n_{mean}\).

Gaussian checks

is_valid_cov(cov[, hbar, rtol, atol])

Checks if the covariance matrix is a valid quantum covariance matrix.

is_pure_cov(cov[, hbar, rtol, atol])

Checks if the covariance matrix is a valid quantum covariance matrix that corresponds to a quantum pure state

is_classical_cov(cov[, hbar, atol])

Checks if the covariance matrix can be efficiently sampled.

fidelity(mu1, cov1, mu2, cov2[, hbar, rtol, ...])

Calculates the fidelity between two Gaussian quantum states.

is_symplectic(S[, rtol, atol])

Checks if the matrix is symplectic.

Conversions

reduced_gaussian(mu, cov, modes)

Returns the vector of means and the covariance matrix of the specified modes.

Xmat(N)

Returns the matrix \(X_n = \begin{bmatrix}0 & I_n\\ I_n & 0\end{bmatrix}\)

Qmat(cov[, hbar])

Returns the \(Q\) Husimi matrix of the Gaussian state.

Covmat(Q[, hbar])

Returns the Wigner covariance matrix in the \(xp\)-ordering of the Gaussian state.

Amat(cov[, hbar, cov_is_qmat])

Returns the \(A\) matrix of the Gaussian state whose hafnian gives the photon number probabilities.

complex_to_real_displacements(mu[, hbar])

Returns the vector of complex displacements and conjugate displacements.

real_to_complex_displacements(beta[, hbar])

Returns the vector of real quadrature displacements.

Means and variances

photon_number_mean(mu, cov, j[, hbar])

Calculate the mean photon number of mode j of a Gaussian state.

photon_number_mean_vector(mu, cov[, hbar])

Calculate the mean photon number of each of the modes in a Gaussian state

photon_number_covar(mu, cov, j, k[, hbar])

Calculate the variance/covariance of the photon number distribution of a Gaussian state.

photon_number_covmat(mu, cov[, hbar])

Calculate the covariance matrix of the photon number distribution of a Gaussian state.

photon_number_expectation(mu, cov, modes[, hbar])

Calculates the expectation value of the product of the number operator of the modes in a Gaussian state.

photon_number_squared_expectation(mu, cov, modes)

Calculates the expectation value of the square of the product of the number operator of the modes in a Gaussian state.

normal_ordered_expectation(mu, cov, rpt[, hbar])

Calculates the expectation value of the normal ordered product \(\prod_{i=0}^{N-1} a_i^{\dagger n_i} \prod_{j=0}^{N-1} a_j^{m_j}\) with respect to an N-mode Gaussian state, where \(\text{rpt}=(n_0, n_1, \ldots, n_{N-1}, m_0, m_1, \ldots, m_{N-1})\).

photon_number_moment(mu, cov, indices[, hbar])

Calculates the expectation value of product of powers of photon number operators of a Gaussian state.

s_ordered_expectation(mu, cov, rpt[, hbar, s])

Calculates the expectation value of the s-ordered product obtained by taking deirvatives of the characteristic function of a Gaussian states, Here, \(\text{rpt}=(n_0, n_1, \ldots, n_{N-1}, m_0, m_1, \ldots, m_{N-1})\).

mean_clicks(cov[, hbar])

Calculates the total mean number of clicks when a zero-mean gaussian state is measured using threshold detectors.

variance_clicks(cov[, hbar])

Calculates the variance of the total number of clicks when a zero-mean gaussian state is measured using threshold detectors.

photon_number_cumulant(mu, cov, modes[, hbar])

Calculates the photon-number cumulant of the modes in the Gaussian state.

click_cumulant(mu, cov, modes[, hbar])

Calculates the click cumulant of the modes in the Gaussian state.

Photon number distributions

pure_state_distribution(cov[, cutoff, hbar, ...])

Calculates the total photon number distribution of a pure state with zero mean.

total_photon_number_distribution(n, k, s, eta)

Probability of observing a total of \(n\) photons when \(k\) identical single-mode squeezed vacua with squeezing parameter \(s\) undergo loss by transmission \(\eta\).

characteristic_function(k, s, eta, mu[, ...])

Calculates the expectation value of the characteristic function \(\langle n^m \exp(mu n) \rangle\) where \(n\) is the total photon number of \(k\) identical single-mode squeezed vacua with squeezing parameter \(s\) undergoing loss by transmission \(\eta\).

Entanglement

vonneumann_entropy(cov[, hbar])

Returns the vonNeumann entropy of a covariance matrix.

entanglement_entropy(cov[, modes_A, split, hbar])

Returns the entanglement entropy of a covariance matrix under a given bipartition.

log_negativity(cov[, modes_A, split, hbar])

Returns the logarithmic negativity of a covariance matrix under a given bipartition.

Code details

Amat(cov, hbar=2, cov_is_qmat=False)[source]

Returns the \(A\) matrix of the Gaussian state whose hafnian gives the photon number probabilities.

Parameters:

  • cov (array) – \(2N\times 2N\) covariance matrix

  • hbar (float) – the value of \(\hbar\) in the commutation relation \([\x,\p]=i\hbar\).

  • cov_is_qmat (bool) – if True, it is assumed that cov is in fact the Q matrix.

Returns:

the \(A\) matrix.

Return type:

array

Covmat(Q, hbar=2)[source]

Returns the Wigner covariance matrix in the \(xp\)-ordering of the Gaussian state. This is the inverse function of Qmat.

Parameters:

  • Q (array) – \(2N\times 2N\) Husimi Q matrix

  • hbar (float) – the value of \(\hbar\) in the commutation relation \([\x,\p]=i\hbar\).

Returns:

the \(xp\)-ordered covariance matrix in the xp-ordering.

Return type:

array

Qmat(cov, hbar=2)[source]

Returns the \(Q\) Husimi matrix of the Gaussian state.

Parameters:

  • cov (array) – \(2N\times 2N xp-\) Wigner covariance matrix

  • hbar (float) – the value of \(\hbar\) in the commutation relation \([\x,\p]=i\hbar\).

Returns:

the \(Q\) matrix.

Return type:

array

Xmat(N)[source]

Returns the matrix \(X_n = \begin{bmatrix}0 & I_n\\ I_n & 0\end{bmatrix}\)

Parameters:

N (int) – positive integer

Returns:

\(2N\times 2N\) array

Return type:

array

adj_scaling(A, n_mean)[source]

Returns the scaling parameter by which the adjacency matrix A should be rescaled so that the Gaussian state that endodes it has a total mean photon number n_mean.

Parameters:

  • A (array) – Adjacency matrix

  • n_mean (float) – Mean photon number of the Gaussian state

Returns:

Scaling parameter

Return type:

float

adj_scaling_torontonian(A, c_mean)[source]

Returns the scaling parameter by which the adjacency matrix A should be rescaled so that the Gaussian state that encodes it has give a mean number of clicks equal to c_mean when measured with threshold detectors.

Parameters:

  • A (array) – adjacency matrix

  • c_mean (float) – mean photon number of the Gaussian state

Returns:

scaling parameter

Return type:

float

adj_to_qmat(A, n_mean)[source]

Returns the Qmat xp-covariance matrix associated to a graph with adjacency matrix \(A\) and with mean photon number \(n_{mean}\).

Parameters:

  • A (array) – a \(N\times N\) np.float64 (symmetric) adjacency matrix

  • n_mean (float) – mean photon number of the Gaussian state

Returns:

the \(2N\times 2N\) Q matrix.

Return type:

array

characteristic_function(k, s, eta, mu, max_iter=10000, delta=1e-14, poly_corr=None)[source]

Calculates the expectation value of the characteristic function \(\langle n^m \exp(mu n) \rangle\) where \(n\) is the total photon number of \(k\) identical single-mode squeezed vacua with squeezing parameter \(s\) undergoing loss by transmission \(\eta\).

Parameters:

  • k (int) – number of squeezed modes

  • s (float) – squeezing parameter

  • eta (float) – transmission parameter, between 0 and 1 inclusive

  • mu (float) – value at which to evaluate the characteristic function

  • max_iter (int) – maximum number of terms allowed in the sum

  • delta (float) – fractional change in the sum after which the sum is stopped

  • poly_corr (int) – give the value of the exponent \(m\) of the polynomial correction

Returns:

the expected value of the moment generation function

Return type:

(float)

complex_to_real_displacements(mu, hbar=2)[source]

Returns the vector of complex displacements and conjugate displacements.

Parameters:

  • mu (array) – length-\(2N\) means vector

  • hbar (float) – the value of \(\hbar\) in the commutation relation \([\x,\p]=i\hbar\).

Returns:

the expectation values \([\langle a_1\rangle, \langle a_2\rangle,\dots,\langle a_N\rangle, \langle a^\dagger_1\rangle, \dots, \langle a^\dagger_N\rangle]\)

Return type:

array

entanglement_entropy(cov, modes_A=None, split=None, hbar=2)[source]

Returns the entanglement entropy of a covariance matrix under a given bipartition.

Parameters:

  • cov (array) – a covariance matrix

  • modes_A (iterable or int) – the subset of modes used for the bipartition

  • split (int) – the index of the mode separating the two partitions (alternative to modes_A)

  • hbar (float) – the value of \(\hbar\) in the commutation relation \([\x,\p]=i\hbar\)

Returns:

logarithmic negativity

Return type:

(float)

fidelity(mu1, cov1, mu2, cov2, hbar=2, rtol=1e-05, atol=1e-08)[source]

Calculates the fidelity between two Gaussian quantum states. For two pure states \(|\psi_1 \rangle, \ |\psi_2 \rangle\) the fidelity is given by \(|\langle \psi_1|\psi_2 \rangle|^2\)

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

The actual implementation used here corresponds to the square of Eq. 96 of ‘Gaussian states and operations - a quick reference’, Brask.

Parameters:

  • mu1 (array) – vector of means of the first state

  • cov1 (array) – covariance matrix of the first state

  • mu2 (array) – vector of means of the second 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:

value of the fidelity between the two states

Return type:

(float)

find_classical_subsystem(cov, hbar=2, atol=1e-08)[source]

Find the largest integer k so that subsystem in modes [0,1,...,k-1] is a classical state.

Parameters:

  • cov (array) – a covariance matrix

  • hbar (float) – value of hbar in the uncertainty relation

  • atol (float) – the absolute tolerance parameter used when determining if the state is classical

Returns:

the largest k so that modes [0,1,...,k-1] are in a classical state.

Return type:

int

fock_tensor(S, alpha, cutoff, choi_r=0.881373587019543, check_symplectic=True, sf_order=False, rtol=1e-05, atol=1e-08)[source]

Calculates the Fock representation of a Gaussian unitary parametrized by the symplectic matrix S and the displacements alpha up to cutoff in Fock space.

Parameters:

  • S (array) – symplectic matrix

  • alpha (array) – complex vector of displacements

  • cutoff (int) – cutoff in Fock space

  • choi_r (float) – squeezing parameter used for the Choi expansion

  • check_symplectic (boolean) – checks whether the input matrix is symplectic

  • sf_order (boolean) – reshapes the tensor so that it follows the sf ordering of indices

  • rtol (float) – the relative tolerance parameter used in np.allclose

  • atol (float) – the absolute tolerance parameter used in np.allclose

Returns:

Tensor containing the Fock representation of the Gaussian unitary

Return type:

(array)

is_classical_cov(cov, hbar=2, atol=1e-08)[source]

Checks if the covariance matrix can be efficiently sampled.

Parameters:

  • cov (array) – a covariance matrix

  • hbar (float) – value of hbar in the uncertainty relation

  • atol (float) – the absolute tolerance parameter used in np.allclose

Returns:

whether the given covariance matrix corresponds to a classical state

Return type:

(bool)

is_pure_cov(cov, hbar=2, rtol=1e-05, atol=1e-08)[source]

Checks if the covariance matrix is a valid quantum covariance matrix that corresponds to a quantum pure state

Parameters:

  • cov (array) – a covariance matrix

  • 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:

whether the given covariance matrix corresponds to a pure state

Return type:

(bool)

is_symplectic(S, rtol=1e-05, atol=1e-08)[source]

Checks if the matrix is symplectic.

Parameters:

  • S (array) – a matrix

  • rtol (float) – the relative tolerance parameter used in np.allclose

  • atol (float) – the absolute tolerance parameter used in np.allclose

Returns:

whether the given matrix is symplectic

Return type:

bool

is_valid_cov(cov, hbar=2, rtol=1e-05, atol=1e-08)[source]

Checks if the covariance matrix is a valid quantum covariance matrix.

Parameters:

  • cov (array) – a covariance matrix

  • 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:

whether the given covariance matrix is a valid covariance matrix

Return type:

(bool)

log_negativity(cov, modes_A=None, split=None, hbar=2)[source]

Returns the logarithmic negativity of a covariance matrix under a given bipartition.

Parameters:

  • cov (array) – a covariance matrix

  • modes_A (iterable or int) – the subset of modes used for the bipartition

  • split (int) – the index of the mode separating the two partitions (alternative to modes_A)

Returns:

entanglement entropy

Return type:

(float)

loss_mat(eta, cutoff)[source]

Constructs a binomial loss matrix with transmission eta up to n photons.

Parameters:

  • eta (float) – Transmission coefficient. eta=0.0 corresponds to complete loss and eta=1.0 corresponds to no loss.

  • cutoff (int) – cutoff in Fock space.

Returns:

\(n\times n\) matrix representing the loss.

Return type:

array

mean_clicks(cov, hbar=2)[source]

Calculates the total mean number of clicks when a zero-mean gaussian state is measured using threshold detectors.

Args

cov (array): \(2N\times 2N\) covariance matrix in xp-ordering hbar (float): the value of \(\hbar\) in the commutation relation \([\x,\p]=i\hbar\)

Returns

float: mean number of clicks

n_body_marginals(mean, cov, cutoff, n, hbar=2)[source]

Calculates the first n-body marginals of a Gaussian state.

For an M-mode Gaussian state there exists a photon number distribution with probability mass function \(p[i_0,i_1,\ldots, i_{M-1}]\). The function n_body_marginals calculates the first n-body marginals of the (all-mode) probability distribution \(p\). The \(n=1\) marginals or single body marginals are simply the probability that mode \(k\) has \(i\) photons, i.e. \(p_k[i]\). For \(n=2\) one obtains the two-body probabilities. For two modes \(k\) and \(l\) this is a two dimensional probability distribution \(p_{k,l}[i,j]\). Note that these marginals have interesting permutation properties, for example \(p_{k,l}[i,j] = p_{l,k}[j,i]\).

The function provided here takes advantage of these symmetries to minimize the amount of calculations. The return of this function is a list of tensors where the first entry contains the one-body marginals of the \(M\) modes (giving a tensor of shape (M, cutoff)), the second entry contains the two-body marginals (giving a tensor of shape (M,M,cutoff, cutoff)) and so on and so forth.

To be clever about not calculating things that can be obtained by permutations it checks whether the index vector representing the modes is sorted. From the way itertools.product works we know that it will always produce a sorted index vector before generating any of its unordered permutations. Thus whenever the index vector is ordered we perform the numerical calculation.

If it is an unsorted index vector it realizes, in the if statement, that it can be obtained by permuting the marginal distribution of something that has already been calculated.

Parameters:

  • mean (array) – length-\(2N\) quadrature displacement vector

  • cov (array) – length-\(2N\) covariance matrix

  • cutoff (int) – cutoff in Fock space

  • n (int) – order of the correlations

  • hbar (float) – the value of \(\hbar\) in the commutation relation \([\x,\p]=i\hbar\).

Returns:

List with arrays containing the \(1,..,n\) body marginal

distributions of the modes

Return type:

list[array]

normal_ordered_expectation(mu, cov, rpt, hbar=2)[source]

Calculates the expectation value of the normal ordered product \(\prod_{i=0}^{N-1} a_i^{\dagger n_i} \prod_{j=0}^{N-1} a_j^{m_j}\) with respect to an N-mode Gaussian state, where \(\text{rpt}=(n_0, n_1, \ldots, n_{N-1}, m_0, m_1, \ldots, m_{N-1})\).

Parameters:

  • mu (array) – length-\(2N\) means vector in xp-ordering.

  • cov (array) – \(2N\times 2N\) covariance matrix in xp-ordering.

  • rpt (list) – integers specifying the terms to calculate.

  • hbar (float) – value of hbar in the uncertainty relation.

Returns:

expectation value of the normal ordered product of operators

Return type:

(float)

photon_number_covar(mu, cov, j, k, hbar=2)[source]

Calculate the variance/covariance of the photon number distribution of a Gaussian state.

Implements the covariance matrix of the photon number distribution of a Gaussian state according to the Last two eq. of Part II. in ‘Multidimensional Hermite polynomials and photon distribution for polymode mixed light’, Dodonov et al.

\[\begin{split}\sigma_{n_j n_j} &= \frac{1}{2}\left(T_j^2 - 2d_j - \frac{1}{2}\right) + \left<\mathbf{Q}_j\right>\mathcal{M}_j\left<\mathbf{Q}_j\right>, \\ \sigma_{n_j n_k} &= \frac{1}{2}\mathrm{Tr}\left(\Lambda_j \mathbf{M} \Lambda_k \mathbf{M}\right) + \left<\mathbf{Q}\right>\Lambda_j \mathbf{M} \Lambda_k\left<\mathbf{Q}\right>,\end{split}\]

where \(T_j\) and \(d_j\) are the trace and the determinant of \(2 \times 2\) matrix \(\mathcal{M}_j\) whose elements coincide with the nonzero elements of matrix \(\mathbf{M}_j = \Lambda_j \mathbf{M} \Lambda_k\) while the two-vector \(\mathbf{Q}_j\) has the components \((q_j, p_j)\). \(2N \times 2N\) projector matrix \(\Lambda_j\) has only two nonzero elements: \(\left(\Lambda_j\right)_{jj} = \left(\Lambda_j\right)_{j+N,j+N} = 1\). Note that the convention for mu used here differs from the one used in Dodonov et al., They both provide the same results in this particular case. Also note that the original reference of Dodonov et al. has an incorrect prefactor of 1/2 in the last terms of the last equation above.

Parameters:

  • mu (array) – vector of means of the Gaussian state using the ordering \([q_1, q_2, \dots, q_n, p_1, p_2, \dots, p_n]\)

  • cov (array) – the covariance matrix of the Gaussian state

  • j (int) – the j th mode

  • k (int) – the k th mode

  • hbar (float) – the hbar convention used in the commutation relation \([q, p]=i\hbar\)

Returns:

the covariance for the photon numbers at modes \(j\) and \(k\).

Return type:

float

photon_number_covmat(mu, cov, hbar=2)[source]

Calculate the covariance matrix of the photon number distribution of a Gaussian state.

Parameters:

  • mu (array) – vector of means of the Gaussian state using the ordering \([q_1, q_2, \dots, q_n, p_1, p_2, \dots, p_n]\)

  • cov (array) – the covariance matrix of the Gaussian state

  • hbar (float) – the hbar convention used in the commutation relation \([q, p]=i\hbar\)

Returns:

the covariance matrix of the photon number distribution

Return type:

array

photon_number_cumulant(mu, cov, modes, hbar=2)[source]

Calculates the photon-number cumulant of the modes in the Gaussian state.

Parameters:

  • mu (array) – length-\(2N\) means vector in xp-ordering.

  • cov (array) – \(2N\times 2N\) covariance matrix in xp-ordering.

  • modes (list or array) – list of modes. Note that it can have repetitions.

  • hbar (float) – value of hbar in the uncertainty relation.

Returns:

the cumulant

Return type:

(float)

photon_number_expectation(mu, cov, modes, hbar=2)[source]

Calculates the expectation value of the product of the number operator of the modes in a Gaussian state.

Parameters:

  • mu (array) – length-\(2N\) means vector in xp-ordering.

  • cov (array) – \(2N\times 2N\) covariance matrix in xp-ordering.

  • modes (list) – list of modes

  • hbar (float) – value of hbar in the uncertainty relation.

Returns:

expectation value of the product of the number operators of the modes.

Return type:

(float)

photon_number_mean(mu, cov, j, hbar=2)[source]

Calculate the mean photon number of mode j of a Gaussian state.

Parameters:

  • mu (array) – vector of means of the Gaussian state using the ordering \([q_1, q_2, \dots, q_n, p_1, p_2, \dots, p_n]\)

  • cov (array) – the covariance matrix of the Gaussian state

  • j (int) – the j th mode

  • hbar (float) – the hbar convention used in the commutation relation \([q, p]=i\hbar\)

Returns:

the mean photon number in mode \(j\).

Return type:

float

photon_number_mean_vector(mu, cov, hbar=2)[source]

Calculate the mean photon number of each of the modes in a Gaussian state

Parameters:

  • mu (array) – vector of means of the Gaussian state using the ordering \([q_1, q_2, \dots, q_n, p_1, p_2, \dots, p_n]\)

  • cov (array) – the covariance matrix of the Gaussian state

  • hbar (float) – the hbar convention used in the commutation relation \([q, p]=i\hbar\)

Returns:

the vector of means of the photon number distribution

Return type:

array

photon_number_moment(mu, cov, indices, hbar=2)[source]

Calculates the expectation value of product of powers of photon number operators of a Gaussian state. The powers are specified by a dictionary with modes as keys and powers as values.

The calculation is performed by first writing any power of the photon number as

\((a^\dagger a)^m = \sum_{k=1}^m c_k a^{\dagger k} a^k\)

where the coefficients \(c_i\) are provided by the function _coeff_normal_ordered.

Parameters:

  • mu (array) – length-\(2N\) means vector in xp-ordering.

  • cov (array) – \(2N\times 2N\) covariance matrix in xp-ordering.

  • indices (dictionary) – specification of the different modes and their power of their photon number

  • hbar (float) – value of hbar in the uncertainty relation.

Returns:

the expectation value of the photon number powers.

Return type:

float

photon_number_squared_expectation(mu, cov, modes, hbar=2)[source]

Calculates the expectation value of the square of the product of the number operator of the modes in a Gaussian state.

Parameters:

  • mu (array) – length-\(2N\) means vector in xp-ordering.

  • cov (array) – \(2N\times 2N\) covariance matrix in xp-ordering.

  • modes (list) – list of modes

  • hbar (float) – value of hbar in the uncertainty relation.

Returns:

expectation value of the square of the product of the number operator of the modes.

Return type:

(float)

probabilities(mu, cov, cutoff, parallel=False, hbar=2.0, rtol=1e-05, atol=1e-08)[source]

Generate the Fock space probabilities of a Gaussian state up to a Fock space cutoff.

Parameters:

  • mu (array) – vector of means of length 2*n_modes

  • cov (array) – covariance matrix of shape [2*n_modes, 2*n_modes]

  • cutoff (int) – cutoff in Fock space

  • parallel (bool) – if True, uses dask for parallelization instead of OpenMP

  • hbar (float) – value of \(\hbar\) in the commutation relation :math;`[hat{x}, hat{p}]=ihbar`

  • rtol (float) – the relative tolerance parameter used in np.allclose

  • atol (float) – the absolute tolerance parameter used in np.allclose

Returns:

Fock space probabilities up to cutoff. The shape of this tensor is [cutoff]*num_modes.

Return type:

(array)

pure_state_distribution(cov, cutoff=50, hbar=2, padding_factor=2)[source]

Calculates the total photon number distribution of a pure state with zero mean.

Parameters:

  • cov (array) – \(2N\times 2N\) covariance matrix in xp-ordering

  • cutoff (int) – Fock cutoff

  • tol (float) – tolerance for determining if displacement is negligible

  • hbar (float) – the value of \(\hbar\) in the commutation

  • padding_factor (int) – expanded size of the photon distribution to avoid accumulation of errors

Returns:

Total photon number distribution

Return type:

(array)

real_to_complex_displacements(beta, hbar=2)[source]

Returns the vector of real quadrature displacements.

Parameters:

  • beta (array) – length-\(2N\) means bivector

  • hbar (float) – the value of \(\hbar\) in the commutation relation \([\x,\p]=i\hbar\).

Returns:

the quadrature expectation values \([\langle q_1\rangle, \langle q_2\rangle,\dots,\langle q_N\rangle, \langle p_1\rangle, \dots, \langle p_N\rangle]\)

Return type:

array

reduced_gaussian(mu, cov, modes)[source]

Returns the vector of means and the covariance matrix of the specified modes.

Parameters:

  • mu (array) – a length-\(2N\) np.float64 vector of means.

  • cov (array) – a \(2N\times 2N\) np.float64 covariance matrix representing an \(N\) mode quantum state.

  • modes (int of Sequence[int]) – indices of the requested modes

Returns:

where means is an array containing the vector of means, and cov is a square array containing the covariance matrix.

Return type:

tuple (means, cov)

s_ordered_expectation(mu, cov, rpt, hbar=2, s=0)[source]

Calculates the expectation value of the s-ordered product obtained by taking deirvatives of the characteristic function of a Gaussian states, Here, \(\text{rpt}=(n_0, n_1, \ldots, n_{N-1}, m_0, m_1, \ldots, m_{N-1})\). indicates how many derivatives are taken with respect to the complex argument and its conjugate. The values \(s=\{1,0,-1\}\) correspond respectively to normal, symmetric and antinormal order.

Parameters:

  • mu (array) – length-\(2N\) means vector in xp-ordering.

  • cov (array) – \(2N\times 2N\) covariance matrix in xp-ordering.

  • rpt (list) – integers specifying the terms to calculate.

  • hbar (float) – value of hbar in the uncertainty relation.

  • s (float) – value setting the ordering it must be between -1 and 1.

Returns:

expectation value of the normal ordered product of operators

Return type:

(float)

total_photon_number_distribution(n, k, s, eta, pref=1.0)[source]

Probability of observing a total of \(n\) photons when \(k\) identical single-mode squeezed vacua with squeezing parameter \(s\) undergo loss by transmission \(\eta\).

For the derivation see Appendix E of ‘Quantum Computational Supremacy via High-Dimensional Gaussian Boson Sampling’, Deshpande et al..

Parameters:

  • n (int) – number of photons

  • k (int) – number of squeezed modes

  • s (float) – squeezing parameter

  • eta (float) – transmission parameter, between 0 and 1 inclusive

  • pref (float) – use to return the probability times pref**n

Returns:

probability of observing a total of n photons or the probability times pref ** n.

Return type:

(float)

tvd_cutoff_bounds(mu, cov, cutoff, hbar=2, check_is_valid_cov=True, rtol=1e-05, atol=1e-08)[source]

Gives bounds of the total variation distance between the exact Gaussian Boson Sampling distribution extending to infinity in Fock space and the ones truncated by any value between 0 and the user provided cutoff.

For the derivation see Appendix B of ‘Exact simulation of Gaussian boson sampling in polynomial space and exponential time’, Quesada and Arrazola.

Parameters:

  • mu (array) – vector of means of the Gaussian state

  • cov (array) – covariance matrix of the Gaussian state

  • cutoff (int) – cutoff in Fock space

  • check_is_valid_cov (bool) – verify that the covariance matrix is physical

  • 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:

values of the bound for different local Fock space dimensions up to cutoff

Return type:

(array)

update_probabilities_with_loss(etas, probs)[source]

Given a list of transmissivities a tensor of probabilitites, calculate an updated tensor of probabilities after loss is applied.

Parameters:

  • etas (list) – List of transmissitivities describing the loss in each of the modes

  • probs (array) – Array of probabilitites in the different modes

Returns:

List of loss-updated probabilities with the same shape as probs.

Return type:

array

update_probabilities_with_noise(probs_noise, probs)[source]

Given a list of noise probability distributions for each of the modes and a tensor of probabilitites, calculate an updated tensor of probabilities after noise is applied.

Parameters:

  • probs_noise (list) – List of probability distributions describing the noise in each of the modes

  • probs (array) – Array of probabilitites in the different modes

Returns:

List of noise-updated probabilities with the same shape as probs.

Return type:

array

variance_clicks(cov, hbar=2)[source]

Calculates the variance of the total number of clicks when a zero-mean gaussian state is measured using threshold detectors.

Args

cov (array): \(2N\times 2N\) covariance matrix in xp-ordering hbar (float): the value of \(\hbar\) in the commutation relation \([\x,\p]=i\hbar\)

Returns

float: variance in the total number of clicks

vonneumann_entropy(cov, hbar=2)[source]

Returns the vonNeumann entropy of a covariance matrix.

Parameters:

  • cov (array) – a covariance matrix

  • hbar (float) – the value of \(\hbar\) in the commutation relation \([\x,\p]=i\hbar\)

Returns:

vonNeumann entropy

Return type:

(float)

Sampling algorithms

Module name: thewalrus.samples

This submodule provides access to algorithms to sample from the hafnian or the torontonian of Gaussian quantum states.

Hafnian sampling

generate_hafnian_sample(cov[, mean, hbar, ...])

Returns a single sample from the Hafnian of a Gaussian state.

hafnian_sample_state(cov, samples[, mean, ...])

Returns samples from the Hafnian of a Gaussian state.

hafnian_sample_graph(A, n_mean[, samples, ...])

Returns samples from the Gaussian state specified by the adjacency matrix \(A\) and with total mean photon number \(n_{mean}\)

hafnian_sample_classical_state(cov, samples)

Returns samples from a Gaussian state that has a positive \(P\) function.

hafnian_sample_graph_rank_one(G, n_mean[, ...])

Returns samples from a rank one adjacency matrix bm{A} = bm{G} bm{G}^T where \(\bm{G}\) is a row vector.

Torontonian sampling

generate_torontonian_sample(cov[, mu, hbar, ...])

Returns a single sample from the Hafnian of a Gaussian state.

torontonian_sample_state(cov, samples[, mu, ...])

Returns samples from the Torontonian of a Gaussian state

torontonian_sample_graph(A, n_mean[, ...])

Returns samples from the Torontonian of a Gaussian state specified by the adjacency matrix \(A\) and with total mean photon number \(n_{mean}\)

torontonian_sample_classical_state(cov, samples)

Returns threshold samples from a Gaussian state that has a positive P function.

threshold_detection_prob(mu, cov, det_pattern)

Threshold detection probabilities for Gaussian states.

Brute force sampling

photon_number_sampler(probabilities, num_samples)

Given a photon-number probability mass function(PMF) it returns samples according to said PMF.

Code details

generate_hafnian_sample(cov, mean=None, hbar=2, cutoff=12, max_photons=8)[source]

Returns a single sample from the Hafnian of a Gaussian state. Code contributed by Jake F.F. Bulmer based on arXiv:2108.01622.

Parameters:

  • cov (array) – a \(2N\times 2N\) np.float64 covariance matrix representing an \(N\) mode quantum state. This can be obtained via the scovmavxp method of the Gaussian backend of Strawberry Fields.

  • mean (array) – a \(2N\) np.float64 vector of means representing the Gaussian state.

  • hbar (float) – (default 2) the value of \(\hbar\) in the commutation relation \([\x,\p]=i\hbar\).

  • cutoff (int) – the Fock basis truncation.

  • max_photons (int) – specifies the maximum number of photons that can be counted.

Returns:

a photon number sample from the Gaussian states.

Return type:

np.array[int]

generate_torontonian_sample(cov, mu=None, hbar=2, max_photons=30, fanout=10, cutoff=1)[source]

Returns a single sample from the Hafnian of a Gaussian state. Code contributed by Jake F.F. Bulmer based on arXiv:2108.01622.

Parameters:

  • cov (array) – a \(2N\times 2N\) np.float64 covariance matrix representing an \(N\) mode quantum state. This can be obtained via the scovmavxp method of the Gaussian backend of Strawberry Fields.

  • mu (array) – a \(2N\) np.float64 displacement vector representing an \(N\) mode quantum state. This can be obtained via the smeanxp method of the Gaussian backend of Strawberry Fields.

  • hbar (float) – (default 2) the value of \(\hbar\) in the commutation relation \([\x,\p]=i\hbar\).

  • max_photons (int) – specifies the maximum number of clicks that can be counted.

Returns:

a threshold sample from the Gaussian state.

Return type:

np.array[int]

hafnian_sample_classical_state(cov, samples, mean=None, hbar=2, atol=1e-08, cutoff=None)[source]

Returns samples from a Gaussian state that has a positive \(P\) function.

Parameters:

  • cov (array) – a \(2N\times 2N\) np.float64 covariance matrix representing an \(N\) mode quantum state. This can be obtained via the scovmavxp method of the Gaussian backend of Strawberry Fields.

  • samples (int) – number of samples to generate

  • mean (array) – vector of means of the gaussian state

  • hbar (float) – the value of \(\hbar\) in the commutation relation \([\x,\p]=i\hbar\).

  • sigdigits (integer) – precision to check that the covariance matrix is a true covariance matrix of a gaussian state.

Returns:

photon number samples from the Gaussian state with covariance cov and vector means mean.

Return type:

np.array[int]

hafnian_sample_graph(A, n_mean, samples=1, cutoff=5, max_photons=30, parallel=False)[source]

Returns samples from the Gaussian state specified by the adjacency matrix \(A\) and with total mean photon number \(n_{mean}\)

Parameters:

  • A (array) – a \(N\times N\) np.float64 (symmetric) adjacency matrix matrix

  • n_mean (float) – mean photon number of the Gaussian state

  • samples (int) – the number of samples to return.

  • cutoff (int) – the Fock basis truncation.

  • max_photons (int) – specifies the maximum number of photons that can be counted.

  • parallel (bool) – if True, uses dask for parallelization of samples

Returns:

photon number samples from the Gaussian state

Return type:

np.array[int]

hafnian_sample_graph_rank_one(G, n_mean, samples=1)[source]

Returns samples from a rank one adjacency matrix bm{A} = bm{G} bm{G}^T where \(\bm{G}\) is a row vector.

Parameters:

  • G (array) – factorization of the rank-one matrix A = G @ G.T.

  • nmean (float) – Total mean photon number.

  • samples (int) – the number of samples to return.

Returns

(array): samples.

hafnian_sample_state(cov, samples, mean=None, hbar=2, cutoff=5, max_photons=30, parallel=False)[source]

Returns samples from the Hafnian of a Gaussian state.

Parameters:

  • cov (array) – a \(2N\times 2N\) np.float64 covariance matrix representing an \(N\) mode quantum state. This can be obtained via the scovmavxp method of the Gaussian backend of Strawberry Fields.

  • samples (int) – the number of samples to return.

  • mean (array) – a \(2N\) np.float64 vector of means representing the Gaussian state.

  • hbar (float) – (default 2) the value of \(\hbar\) in the commutation relation \([\x,\p]=i\hbar\).

  • cutoff (int) – the Fock basis truncation.

  • max_photons (int) – specifies the maximum number of photons that can be counted.

  • parallel (bool) – if True, uses dask for parallelization of samples

Returns:

photon number samples from the Gaussian state

Return type:

np.array[int]

photon_number_sampler(probabilities, num_samples, out_of_bounds=False)[source]

Given a photon-number probability mass function(PMF) it returns samples according to said PMF.

Parameters:

  • probabilities (array) – probability tensor of the modes, has shape [cutoff]*num_modes

  • num_samples (int) – number of samples requested

  • out_of_bounds (boolean) – if False the probability distribution is renormalized. If not False, the value of out_of_bounds is used as a placeholder for samples where more than the cutoff of probabilities are detected.

Returns:

Samples, with shape [num_sample, num_modes]

Return type:

(array)

threshold_detection_prob(mu, cov, det_pattern, hbar=2, atol=1e-10, rtol=1e-10)[source]

Threshold detection probabilities for Gaussian states. Formula from Jake Bulmer, Nicolas Quesada and Stefano Paesani. When state is displaced, threshold_detection_prob_displacement is called. Otherwise, tor is called.

Parameters:

  • mu (1d array) – means of xp Gaussian Wigner function

  • cov (2d array) – : xp Wigner covariance matrix

  • det_pattern (1d array) – array of {0,1} to describe the threshold detection outcome

  • hbar (float) – the value of \(\hbar\) in the commutation relation \([\x,\p]=i\hbar\).

  • rtol (float) – the relative tolerance parameter used in np.allclose

  • atol (float) – the absolute tolerance parameter used in np.allclose

Returns:

probability of detection pattern

Return type:

np.float64

torontonian_sample_classical_state(cov, samples, mean=None, hbar=2, atol=1e-08)[source]

Returns threshold samples from a Gaussian state that has a positive P function.

Parameters:

  • cov (array) – a \(2N\times 2N\) np.float64 covariance matrix representing an \(N\) mode quantum state. This can be obtained via the scovmavxp method of the Gaussian backend of Strawberry Fields.

  • samples (int) – number of samples to generate

  • mean (array) – vector of means of the Gaussian state

  • hbar (float) – the value of \(\hbar\) in the commutation relation \([\x,\p]=i\hbar\).

  • sigdigits (integer) – precision to check that the covariance matrix is a true covariance matrix of a gaussian state.

Returns:

threshold samples from the Gaussian state with covariance cov and vector means mean.

Return type:

np.array[int]

torontonian_sample_graph(A, n_mean, samples=1, max_photons=30, fanout=10, cutoff=1, parallel=False)[source]

Returns samples from the Torontonian of a Gaussian state specified by the adjacency matrix \(A\) and with total mean photon number \(n_{mean}\)

Parameters:

  • A (array) – a \(N\times N\) np.float64 (symmetric) adjacency matrix matrix

  • n_mean (float) – mean photon number of the Gaussian state

  • samples (int) – the number of samples to return.

  • max_photons (int) – specifies the maximum number of photons that can be counted.

  • parallel (bool) – if True, uses dask for parallelization of samples

Returns:

photon number samples from the Torontonian of the Gaussian state

Return type:

np.array[int]

torontonian_sample_state(cov, samples, mu=None, hbar=2, max_photons=30, fanout=10, cutoff=1, parallel=False)[source]

Returns samples from the Torontonian of a Gaussian state

Parameters:

  • cov (array) – a \(2N\times 2N\) np.float64 covariance matrix representing an \(N\) mode quantum state. This can be obtained via the scovmavxp method of the Gaussian backend of Strawberry Fields.

  • samples (int) – number of samples to generate

  • mu (array) – a \(2N\) np.float64 displacement vector representing an \(N\) mode quantum state. This can be obtained via the smeanxp method of the Gaussian backend of Strawberry Fields.

  • hbar (float) – (default 2) the value of \(\hbar\) in the commutation relation \([\x,\p]=i\hbar\).

  • max_photons (int) – specifies the maximum number of clicks that can be counted.

  • parallel (bool) – if True, uses dask for parallelization of samples

Returns:

threshold samples from the Gaussian state.

Return type:

np.array[int]

Classical sampling algorithms

Module name: thewalrus.csamples

This submodule provides access to classical sampling algorithms for thermal states going through an interferometer specified by a real orthogonal matrix. The quantum state to be sampled is specified by a positive semidefinite real matrix and a mean photon number. The algorithm implemented here was first derived in

For more precise details of the implementation see

  • Soran Jahangiri, Juan Miguel Arrazola, Nicolás Quesada, and Nathan Killoran. “Point processes with Gaussian boson sampling” Phys. Rev. E 101, 022134, (2020)..

Summary

rescale_adjacency_matrix_thermal(A, n_mean)

Returns the scaling parameter by which the adjacency matrix A should be rescaled so that the Gaussian state that encodes it has a total mean photon number n_mean for thermal sampling.

rescale_adjacency_matrix(A, n_mean, scale[, ...])

Returns the scaling parameter by which the adjacency matrix A should be rescaled so that the Gaussian state that encodes it has a total mean photon number n_mean.

generate_thermal_samples(ls, O[, num_samples])

Generates samples of the Gaussian state in terms of the mean photon number parameter ls and the interferometer O.

Code details

generate_thermal_samples(ls, O, num_samples=1)[source]

Generates samples of the Gaussian state in terms of the mean photon number parameter ls and the interferometer O.

Parameters:

  • ls (array) – squashing parameters

  • O (array) – Orthogonal matrix representing the interferometer

  • num_samples – Number of samples to generate

Returns:

samples

Return type:

list(array

rescale_adjacency_matrix(A, n_mean, scale, check_positivity=True, check_symmetry=True, rtol=1e-05, atol=1e-08)[source]

Returns the scaling parameter by which the adjacency matrix A should be rescaled so that the Gaussian state that encodes it has a total mean photon number n_mean.

Parameters:

  • A (array) – Adjacency matrix, assumed to be positive semi-definite and real

  • n_mean (float) – Mean photon number of the Gaussian state

  • scale (float) – Determines whether to rescale the matrix for thermal sampling (scale = 1.0) or for squashed sampling (scale = 2.0)

  • check_positivity (bool) – Checks if the matrix A is positive semidefinite

  • check_symmetry (bool) – Checks if the matrix is symmetric

  • rtol – relative tolerance for the checks

  • atol – absolute tolerance for the checks

Returns:

rescaled eigenvalues and eigenvectors of the matrix A

Return type:

tuple(array,array)

rescale_adjacency_matrix_thermal(A, n_mean, check_positivity=True, check_symmetry=True, rtol=1e-05, atol=1e-08)[source]

Returns the scaling parameter by which the adjacency matrix A should be rescaled so that the Gaussian state that encodes it has a total mean photon number n_mean for thermal sampling.

Parameters:

  • A (array) – Adjacency matrix, assumed to be positive semi-definite and real

  • n_mean (float) – Mean photon number of the Gaussian state

  • check_positivity (bool) – Checks if the matrix A is positive semidefinite

  • check_symmetry (bool) – Checks if the matrix is symmetric

  • rtol – relative tolerance for the checks

  • atol – absolute tolerance for the checks

Returns:

rescaled eigenvalues and eigenvectors of the matrix A

Return type:

tuple(array,array)

Symplectic operations

Module name: thewalrus.symplectic

Contains some Gaussian operations and auxiliary functions.

Auxiliary functions

expand(S, modes, N)

Expands a Symplectic matrix S to act on the entire subsystem.

expand_vector(alpha, mode, N[, hbar])

Returns the phase-space displacement vector associated to a displacement.

expand_passive(T, modes, N)

Returns the expanded linear optical transformation acting on specified modes, with identity acting on all other modes

reduced_state(mu, cov, modes)

Returns the vector of means and the covariance matrix of the specified modes.

is_symplectic(S[, rtol, atol])

Checks if matrix S is a symplectic matrix

sympmat(N[, dtype])

Returns the matrix \(\Omega_n = \begin{bmatrix}0 & I_n\\ -I_n & 0\end{bmatrix}\)

xxpp_to_xpxp(S)

Permutes the entries of the input from xxpp ordering to xpxp ordering.

xpxp_to_xxpp(S)

Permutes the entries of the input from xpxp ordering to xxpp ordering.

Gaussian states

vacuum_state(modes[, hbar, dtype])

Returns the vacuum state.

Gates and operations

two_mode_squeezing(r, phi[, dtype])

Two-mode squeezing.

squeezing(r[, phi, dtype])

Squeezing.

interferometer(U)

Interferometer.

loss(mu, cov, T, mode[, nbar, hbar])

Loss channel acting on a Gaussian state.

mean_photon_number(mu, cov[, hbar])

Calculates the mean photon number for a given one-mode state.

beam_splitter(theta, phi[, dtype])

Beam-splitter.

rotation(theta[, dtype])

Rotation gate.

Code details

beam_splitter(theta, phi, dtype=<class 'numpy.float64'>)[source]

Beam-splitter.

Parameters:

  • theta (float) – transmissivity parameter

  • phi (float) – phase parameter

  • dtype (numpy.dtype) – datatype to represent the Symplectic matrix

Returns:

symplectic-orthogonal transformation matrix of an interferometer with angles theta and phi

Return type:

array

expand(S, modes, N)[source]

Expands a Symplectic matrix S to act on the entire subsystem. If the input is a single mode symplectic, then extends it to act on multiple modes.

Supports scipy sparse matrices. Instances of coo_array, dia_array, bsr_array will be transformed into csr_array`.

Parameters:

  • S (ndarray or spmatrix) – a \(2M\times 2M\) Symplectic matrix

  • modes (Sequence[int]) – the list of modes S acts on

  • N (int) – full size of the subsystem

Returns:

the resulting \(2N\times 2N\) Symplectic matrix

Return type:

array

expand_passive(T, modes, N)[source]

Returns the expanded linear optical transformation acting on specified modes, with identity acting on all other modes

Parameters:

  • T (array) – square \(M \times M\) matrix of linear optical transformation

  • modes (array) – the \(M\) modes of the transformation

  • N (int) – number of modes in the new expanded transformation

Returns:

\(N \times N\) array of expanded passive transformation

Return type:

array

expand_vector(alpha, mode, N, hbar=2.0)[source]

Returns the phase-space displacement vector associated to a displacement.

Parameters:

  • alpha (complex) – complex displacement

  • mode (int) – mode index

  • N (int) – number of modes

Returns:

phase-space displacement vector of size 2*N

Return type:

array

interferometer(U)[source]

Interferometer.

Parameters:

U (array) – unitary matrix

Returns:

symplectic transformation matrix

Return type:

array

is_symplectic(S, rtol=1e-05, atol=1e-08)[source]

Checks if matrix S is a symplectic matrix

Parameters:

S (array) – a square matrix

Returns:

whether the given matrix is symplectic

Return type:

(bool)

loss(mu, cov, T, mode, nbar=0, hbar=2)[source]

Loss channel acting on a Gaussian state.

Parameters:

  • mu (array) – means vector

  • cov (array) – covariance matri

  • T (float) – transmission; 1 corresponds to no loss, 0 to full loss.

  • mode (int) – mode to act on

  • nbar (float) – thermal mean population (default 0)

  • hbar (float) – (default 2) the value of \(\hbar\) in the commutation relation \([\x,\p]=i\hbar\)

Returns:

the means vector and covariance matrix of the resulting state

Return type:

tuple[array]

mean_photon_number(mu, cov, hbar=2)[source]

Calculates the mean photon number for a given one-mode state.

Parameters:

  • mu (array) – length-2 vector of means

  • cov (array) – \(2\times 2\) covariance matrix

  • hbar (float) – (default 2) the value of \(\hbar\) in the commutation relation \([\x,\p]=i\hbar\)

Returns:

the photon number expectation and variance

Return type:

tuple

passive_transformation(mu, cov, T, hbar=2)[source]

Perform a covariance matrix transformation for an arbitrary linear optical channel on an \(N\) modes state mapping it to a to an \(M\) modes state.

Parameters:

  • mu (array) – \(2N\)-length means vector

  • cov (array) – \(2N \times 2N\) covariance matrix

  • T (array) – \(M \times N\) linear optical transformation

Keyword Arguments:

hbar (float) – the value to use for hbar

Returns:

\(2M\)-length transformed means vector array \(2M \times 2M\) tranformed covariance matrix

Return type:

array

reduced_state(mu, cov, modes)[source]

Returns the vector of means and the covariance matrix of the specified modes.

Parameters:

modes (int of Sequence[int]) – indices of the requested modes

Returns:

means is an array containing the vector of means, and cov is a square array containing the covariance matrix. Both use the \(xp\)-ordering.

Return type:

tuple (means, cov)

rotation(theta, dtype=<class 'numpy.float64'>)[source]

Rotation gate.

Parameters:

  • theta (float) – rotation angle

  • dtype (numpy.dtype) – datatype to represent the Symplectic matrix

Returns:

rotation matrix by angle theta

Return type:

array

squeezing(r, phi=None, dtype=<class 'numpy.float64'>)[source]

Squeezing. In fock space this corresponds to:

\[\exp(\tfrac{1}{2}r e^{i \phi} (a^2 - a^{\dagger 2}) ).\]

By passing an array of squeezing parameters and phases, it applies a tensor product of squeezing operations.

Parameters:

  • r (Union[array, float]) – squeezing magnitude

  • phi (Union[array, float]) – rotation parameter. If None, then the function uses zeros of the same shape as r.

  • dtype (numpy.dtype) – datatype to represent the Symplectic matrix. Defaults to numpy.float64.

Returns:

symplectic transformation matrix

Return type:

array

sympmat(N, dtype=<class 'numpy.float64'>)[source]

Returns the matrix \(\Omega_n = \begin{bmatrix}0 & I_n\\ -I_n & 0\end{bmatrix}\)

Parameters:

  • N (int) – positive integer

  • dtype (numpy.dtype) – datatype to represent the Symplectic matrix

Returns:

\(2N\times 2N\) array

Return type:

array

two_mode_squeezing(r, phi, dtype=<class 'numpy.float64'>)[source]

Two-mode squeezing.

Parameters:

  • r (float) – squeezing magnitude

  • phi (float) – rotation parameter

  • dtype (numpy.dtype) – datatype to represent the Symplectic matrix

Returns:

symplectic transformation matrix

Return type:

array

vacuum_state(modes, hbar=2.0, dtype=<class 'numpy.float64'>)[source]

Returns the vacuum state.

Parameters:

  • modes (str) – Returns the vector of means and the covariance matrix

  • hbar (float) – (default 2) the value of \(\hbar\) in the commutation relation \([\x,\p]=i\hbar\)

  • dtype (numpy.dtype) – datatype to represent the covariance matrix and vector of means

Returns:

the means vector and covariance matrix of the vacuum state

Return type:

list[array]

xpxp_to_xxpp(S)[source]

Permutes the entries of the input from xpxp ordering to xxpp ordering.

Parameters:

S (array) – input even dimensional square matrix or vector

Returns:

permuted matrix or vector

Return type:

(array)

xxpp_to_xpxp(S)[source]

Permutes the entries of the input from xxpp ordering to xpxp ordering.

Parameters:

S (array) – input even dimensional square matrix or array

Returns:

permuted matrix or array

Return type:

(array)

Characteristic polynomials

Module name: thewalrus.charpoly

This module implements La Budde’s algorithm to calculate the characteristic polynomials of matrices.

Summary

get_reflection_vector(matrix, k)

Compute reflection vector for householder transformation on general complex matrices.

apply_householder(A, v, k)

Apply householder transformation on a matrix A.

reduce_matrix_to_hessenberg(matrix)

Reduce the matrix to upper hessenberg form without Lapack.

charpoly(H)

Calculates the characteristic polynomial of the matrix H.

powertrace(H, n)

Calculates the powertraces of the matrix H up to power n-1.

Code details

apply_householder(A, v, k)[source]

Apply householder transformation on a matrix A. See Matrix Computations by Golub and Van Loan (4th Edition) Sections 5.1.4 and 7.4.2.

Parameters:

  • A (array) – A matrix to apply householder on

  • v (array) – reflection vector

  • size_A (int) – size of matrix A

  • k (int) – offset for submatrix

charpoly(H)[source]

Calculates the characteristic polynomial of the matrix H.

Parameters:

H (array) – square matrix

Returns

(array): list of power traces from 0 to n-1

get_reflection_vector(matrix, k)[source]

Compute reflection vector for householder transformation on general complex matrices. See Introduction to Numerical Analysis-Springer New York (2002) (3rd Edition) by J. Stoer and R. Bulirsch Section 6.5.1.

Parameters:

  • matrix (array) – the matrix in the householder transformation

  • k (int) – offset for submatrix

Returns:

reflection vector

Return type:

array

powertrace(H, n)[source]

Calculates the powertraces of the matrix H up to power n-1.

Parameters:

  • H (array) – square matrix

  • n (int) – required order

Returns:

list of power traces from 0 to n-1

Return type:

(array)

reduce_matrix_to_hessenberg(matrix)[source]

Reduce the matrix to upper hessenberg form without Lapack. This function only accepts Row-Order matrices.

Parameters:

matrix (array) – the matrix to be reduced

Returns:

matrix in hessenberg form

Return type:

array

Random matrices

Module name: thewalrus.random

This submodule provides access to utility functions to generate random unitary, symplectic and covariance matrices.

Summary

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.

random_block_interferometer(N[, top_one, real])

Generates a random interferometer with blocks of at most size 2.

random_banded_interferometer(N, w[, ...])

Generates a banded unitary matrix.

Code details

randnc(*arg)[source]

Normally distributed array of random complex numbers.

random_banded_interferometer(N, w, top_one_init=True, real=False)[source]

Generates a banded unitary matrix.

Parameters:

  • N (int) – number of modes

  • w (int) – bandwidth

  • top_one_init (bool) – if True places a 1times1 interferometer in the top-left-most block of the first matrix in the product

  • real (bool) – return a random real orthogonal matrix

Returns:

random \(N\times N\) unitary with the specified block structure

Return type:

array

random_block_interferometer(N, top_one=True, real=False)[source]

Generates a random interferometer with blocks of at most size 2.

Parameters:

  • N (int) – number of modes

  • top_one (bool) – if True places a 1times1 interferometer in the top-left most block

  • real (bool) – return a random real orthogonal matrix

Returns:

random \(N\times N\) unitary with the specified block structure

Return type:

array

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 \(x\) 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

random_interferometer(N, real=False)[source]

Random unitary matrix representing an interferometer. For more details, see [41].

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

random_symplectic(N, passive=False, block_diag=False, scale=1.0)[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. Note that for the Symplectic group there is no notion of Haar measure since this is group is not compact.

Parameters:

  • N (int) – number of modes

  • passive (bool) – If True, returns a passive Gaussian transformation (i.e., one that preserves photon number). If False, 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

  • scale (float) – Sets the scale of the random values used as squeezing parameters. They will range from 0 to \(\sqrt{2}\texttt{scale}\)

Returns:

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

Return type:

array

Fock representations

Module name: thewalrus.fock_gradients

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

Summary

displacement(r, phi, cutoff[, dtype])

Calculates the matrix elements of the displacement gate using a recurrence relation.

squeezing(r, theta, cutoff[, dtype])

Calculates the matrix elements of the squeezing gate using a recurrence relation.

beamsplitter(theta, phi, cutoff[, dtype])

Calculates the Fock representation of the beamsplitter.

two_mode_squeezing(r, theta, cutoff[, dtype])

Calculates the matrix elements of the two-mode squeezing gate recursively.

mzgate(theta, phi, cutoff[, dtype])

Calculates the Fock representation of the Mach-Zehnder interferometer.

grad_displacement(T, r, phi)

Calculates the gradients of the displacement gate with respect to the displacement magnitude and angle.

grad_squeezing(T, r, phi)

Calculates the gradients of the squeezing gate with respect to the squeezing magnitude and angle

grad_beamsplitter(T, theta, phi)

Calculates the gradients of the beamsplitter gate with respect to the transmissivity angle and reflection phase

grad_two_mode_squeezing(T, r, theta)

Calculates the gradients of the two-mode squeezing gate with respect to the squeezing magnitude and angle

grad_mzgate(T, theta, phi)

Calculates the gradients of the Mach-Zehnder interferometer with respect to the transmissivity angle and reflection phase

Code details

beamsplitter(theta, phi, cutoff, dtype=<class 'numpy.complex128'>)[source]

Calculates the Fock representation of the beamsplitter.

Parameters:

  • theta (float) – transmissivity angle of the beamsplitter. The transmissivity is \(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:

The Fock representation of the gate

Return type:

array[float]

displacement(r, phi, cutoff, dtype=<class 'numpy.complex128'>)[source]

Calculates the matrix elements of the displacement gate using a recurrence relation.

Parameters:

  • 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:

matrix representing the displacement operation.

Return type:

array[complex]

grad_beamsplitter(T, theta, phi)[source]

Calculates the gradients of the beamsplitter gate with respect to the transmissivity angle and reflection phase

Parameters:

  • T (array[complex]) – array representing the gate

  • theta (float) – transmissivity angle of the beamsplitter. The transmissivity is \(t=\cos(\theta)\)

  • phi (float) – reflection phase of the beamsplitter

Returns:

The gradient of the beamsplitter gate with respect to theta and phi

Return type:

tuple[array[complex], array[complex]]

grad_displacement(T, r, phi)[source]

Calculates the gradients of the displacement gate with respect to the displacement magnitude and angle.

Parameters:

  • T (array[complex]) – array representing the gate

  • r (float) – displacement magnitude

  • phi (float) – displacement angle

Returns:

The gradient of the displacement gate with respect to r and phi

Return type:

tuple[array[complex], array[complex]]

grad_mzgate(T, theta, phi)[source]

Calculates the gradients of the Mach-Zehnder interferometer with respect to the transmissivity angle and reflection phase

Parameters:

  • T (array[complex]) – array representing the gate

  • theta (float) – internal of the mzgate

  • phi (float) – external phase of the mzgate

Returns:

The gradient of the mzgate gate with respect to theta and phi

Return type:

tuple[array[complex], array[complex]]

grad_squeezing(T, r, phi)[source]

Calculates the gradients of the squeezing gate with respect to the squeezing magnitude and angle

Parameters:

  • T (array[complex]) – array representing the gate

  • r (float) – squeezing magnitude

  • phi (float) – squeezing angle

Returns:

The gradient of the squeezing gate with respect to the r and phi

Return type:

tuple[array[complex], array[complex]]

grad_two_mode_squeezing(T, r, theta)[source]

Calculates the gradients of the two-mode squeezing gate with respect to the squeezing magnitude and angle

Parameters:

  • T (array[complex]) – array representing the gate

  • r (float) – squeezing magnitude

  • theta (float) – squeezing angle

Returns:

The gradient of the two-mode squeezing gate with respect to r and phi

Return type:

tuple[array[complex], array[complex]]

mzgate(theta, phi, cutoff, dtype=<class 'numpy.complex128'>)[source]

Calculates the Fock representation of the Mach-Zehnder interferometer.

Parameters:

  • 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:

The Fock representation of the gate

Return type:

array[float]

squeezing(r, theta, cutoff, dtype=<class 'numpy.complex128'>)[source]

Calculates the matrix elements of the squeezing gate using a recurrence relation.

Parameters:

  • 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:

matrix representing the squeezing gate.

Return type:

array[complex]

two_mode_squeezing(r, theta, cutoff, dtype=<class 'numpy.complex128'>)[source]

Calculates the matrix elements of the two-mode squeezing gate recursively.

Parameters:

  • 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:

The Fock representation of the gate

Return type:

array[float]

Decompositions

Module name: thewalrus.decompositions

This module implements common shared matrix decompositions that are used to perform gate decompositions.

For mathematical details of these decompositions see

Houde et al. Matrix decompositions in Quantum Optics: Takagi/Autonne, Bloch-Messiah/Euler, Iwasawa, and Williamson [42] Summary ——-

williamson(V[, rtol, atol])

Williamson decomposition of positive-definite (real) symmetric matrix.

symplectic_eigenvals(cov)

Returns the symplectic eigenvalues of a covariance matrix.

blochmessiah(S)

Returns the Bloch-Messiah decomposition of a symplectic matrix S = O @ D @ Q

takagi(A[, svd_order])

Autonne-Takagi decomposition of a complex symmetric (not Hermitian!) matrix.

pre_iwasawa(S)

Pre-Iwasawa decomposition of a symplectic matrix.

iwasawa(S)

Iwasawa decomposition of a symplectic matrix.

Code details

blochmessiah(S)[source]
Returns the Bloch-Messiah decomposition of a symplectic matrix S = O @ D @ Q

where O and Q are orthogonal symplectic matrices and D is a positive-definite diagonal matrix of the form diag(d1,d2,…,dn,d1^-1, d2^-1,…,dn^-1).

See Houde et al. Matrix decompositions in Quantum Optics: Takagi/Autonne, Bloch-Messiah/Euler, Iwasawa, and Williamson

Parameters:

S (array[float]) – 2N x 2N real symplectic matrix

Returns:

orthogonal symplectic matrix O

array[float], : diagonal matrix D array[float]) : orthogonal symplectic matrix Q

Return type:

tuple(array[float],

iwasawa(S)[source]

Iwasawa decomposition of a symplectic matrix. See Arvind et al. The Real Symplectic Groups in Quantum Mechanics and Optics and Houde et al. Matrix decompositions in Quantum Optics: Takagi/Autonne, Bloch-Messiah/Euler, Iwasawa, and Williamson

Parameters:

S (array) – the symplectic matrix

Returns:

(E,D,F) symplectic matrices such that E @ D @ F = S, EE = np.block([[AA, np.zeros(N,N)],[CC, np.linalg.inv(A.T)]]) with A.T @ C == C.T @ A, and AA upper trinagular with ones in the diagonal DD is diagonal and symplectic, FF is symplectic orthogonal.

Return type:

tuple[array, array, array]

pre_iwasawa(S)[source]

Pre-Iwasawa decomposition of a symplectic matrix. See Arvind et al. The Real Symplectic Groups in Quantum Mechanics and Optics and Houde et al. Matrix decompositions in Quantum Optics: Takagi/Autonne, Bloch-Messiah/Euler, Iwasawa, and Williamson

Parameters:

S (array) – the symplectic matrix

Returns:

(E,D,F) symplectic matrices such that E @ D @ F = S and, E = np.block([[np.eye(N), np.zeros(N,N)],[X, np.eye(N)]]) with X == X.T, D is block diagonal with the top left block being the inverse of the bottom right block, F is symplectic orthogonal.

Return type:

tuple[array, array, array]

symplectic_eigenvals(cov)[source]

Returns the symplectic eigenvalues of a covariance matrix.

Parameters:

cov (array) – a covariance matrix

Returns:

symplectic eigenvalues

Return type:

(array)

takagi(A, svd_order=True)[source]

Autonne-Takagi decomposition of a complex symmetric (not Hermitian!) matrix. Note that the input matrix is internally symmetrized by taking its upper triangular part. If the input matrix is indeed symmetric this leaves it unchanged.

See Houde et al. Matrix decompositions in Quantum Optics: Takagi/Autonne, Bloch-Messiah/Euler, Iwasawa, and Williamson

Parameters:

  • A (array) – square, symmetric matrix

  • svd_order (boolean) – whether to return result by ordering the singular values of A in descending (True) or ascending (False) order.

Returns:

(r, U), where r are the singular values, and U is the Autonne-Takagi unitary, such that \(A = U \diag(r) U^T\).

Return type:

tuple[array, array]

williamson(V, rtol=1e-05, atol=1e-08)[source]

Williamson decomposition of positive-definite (real) symmetric matrix.

See Houde et al. Matrix decompositions in Quantum Optics: Takagi/Autonne, Bloch-Messiah/Euler, Iwasawa, and Williamson

Parameters:

  • V (array[float]) – positive definite symmetric (real) matrix

  • rtol (float) – the relative tolerance parameter used in np.allclose

  • atol (float) – the absolute tolerance parameter used in np.allclose

Returns:

(Db, S) where Db is a diagonal matrix

and S is a symplectic matrix such that \(V = S Db S^T\)

Return type:

tuple[array,array]

Reference implementations

Module name: thewalrus.reference

This submodule provides access to reference implementations of the hafnian and loop hafnian by summing over the set of perfect matching permutations or the set of single pair matchings.

For more details on these definitions see:

  • Andreas Björklund, Brajesh Gupt, and Nicolás Quesada. “A faster hafnian formula for complex matrices and its benchmarking on the Titan supercomputer” arxiv:1805.12498 (2018)

Reference functions

hafnian(M[, loop])

Returns the (loop) hafnian of the matrix \(M\).

Code details

Auxiliary functions

memoized(f[, maxsize])

Decorator used to memoize a generator.

partitions(s[, singles, pairs])

Returns the partitions of a tuple in terms of pairs and singles.

spm(s)

Generator for the set of single pair matchings.

pmp(s)

Generator for the set of perfect matching permutations.

rspm(s)

Generates the restricted set of single-pair matchings.

rpmp(s)

Generates the restricted set of perfect matchings matching permutations.

T(n)

Returns the \(n\) th telephone number.

Code details

T(n)[source]

Returns the \(n\) th telephone number.

They satisfy the recursion relation \(T(n) = T(n-1)+(n-1)T(n-2)\) and \(T(0)=T(1)=1\).

See https://oeis.org/A000085 for more details.

Parameters:

n (integer) – index

Returns:

the nth telephone number

Return type:

int

bitstrings(n)[source]

Returns the bistrings from 0 to n/2

Parameters:

n (int) – Twice the highest bitstring value.

Returns:

An iterable of all bistrings.

Return type:

(iterator)

mapper(x, objects)[source]

Helper function to turn a permutation and bistring into an element of rpmp.

Parameters:

  • x (tuple) – tuple containing a permutation and a bistring

  • objects (list) – list objects to permute

Returns:

permuted objects

Return type:

tuple

memoized(f, maxsize=1000)[source]

Decorator used to memoize a generator.

The standard approach of using functools.lru_cache cannot be used, as it only memoizes the generator object, not the results of the generator.

See https://stackoverflow.com/a/10726355/10958457 for the original implementation.

Parameters:

  • f (function or generator) – function or generator to memoize

  • maxsize (int) – positive integer that defines the maximum size of the cache

Returns:

the memoized function or generator

Return type:

function or generator

mtl(A, loop=False)[source]

Returns the Montrealer of an NxN matrix and an N-length vector.

Parameters:

  • A (array) – an NxN array of even dimensions. Can be symbolic.

  • loop (boolean) – if set to True, the loop montrealer is returned

Returns:

the Montrealer of matrix A.

Return type:

np.float64, np.complex128 or sympy.core.add.Add

partitions(s, singles=True, pairs=True)[source]

Returns the partitions of a tuple in terms of pairs and singles.

Parameters:

  • s (tuple) – a tuple representing the (multi-)set that will be partitioned. Note that it must hold that len(s) >= 3.

  • single (boolean) – allow singles in the partitions

  • pairs (boolean) – allow pairs in the partitions

Returns:

a generators that goes through all the single-double partitions of the tuple

Return type:

generator

pmp(s)[source]

Generator for the set of perfect matching permutations.

Parameters:

s (tuple) – an input tuple

Returns:

the set of perfect matching permutations of the tuple s

Return type:

generator

rpmp(s)[source]

Generates the restricted set of perfect matchings matching permutations.

Parameters:

s (tuple) – tuple of labels to be used

Returns:

the set of restricted perfect matching permutations of the tuple s

Return type:

generator

rspm(s)[source]

Generates the restricted set of single-pair matchings.

Parameters:

s (tuple) – tuple of labels to be used

Returns:

the set of restricted perfect matching permutations of the tuple s

Return type:

generator

splitter(elem)[source]

Takes an element from the restricted perfect matching permutations (rpmp) and returns all the associated elements in the restricted single pair matchings (rspm)

Parameters:

elem (tuple) – tuple representing an element of rpmp

Returns:

all the associated elements in rspm

Return type:

(iterator)

spm(s)[source]

Generator for the set of single pair matchings.

Parameters:

s (tuple) – an input tuple

Returns:

the set of single pair matching of the tuple s

Return type:

generator

Contents

Getting started

Background

The Walrus API