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