Python library

This is the top level module of the The Walrus Python interface, 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, [6]. 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 [5]. This algorithm scales like \(\mathcal{O}(n^4 2^{n/2})\). This algorithm does not currently support the loop hafnian.

Repeating hafnian algorithm

The algorithm described in From moments of sum to moments of product, [2]. This method is more efficient for matrices with repeated rows and columns, and supports caclulation 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, [14]. 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 [18].

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, [6].

Python wrappers

hafnian(A[, loop, recursive, tol, quad, …])

Returns the hafnian of a matrix.

hafnian_repeated(A, rpt[, mu, loop, tol])

Returns the hafnian of matrix with repeated rows/columns.

hafnian_batched(A, cutoff[, mu, tol, …])

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

tor(A[, fsum])

Returns the Torontonian of a matrix.

perm(A[, quad, fsum])

Returns the permanent of a matrix via the Ryser formula.

permanent_repeated(A, rpt)

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

hermite_multidimensional(R, cutoff[, y, …])

Returns the multidimensional Hermite polynomials \(H_k^{(R)}(y)\).

Pure Python functions

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

hafnian(A, loop=False, recursive=True, tol=1e-12, quad=True, approx=False, num_samples=1000)[source]

Returns the hafnian of a matrix.

For more direct control, you may wish to call haf_real(), haf_complex(), or haf_int() directly.

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

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

  • recursive (bool) – If True, the recursive algorithm is used. Note: the recursive algorithm does not currently support the loop hafnian. If loop=True, then this keyword argument is ignored.

  • tol (float) – the tolerance when checking that the matrix is symmetric. Default tolerance is 1e-12.

  • quad (bool) – If True, the hafnian algorithm is performed with quadruple precision.

  • 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

np.int64 or np.float64 or np.complex128

hafnian_repeated(A, rpt, mu=None, loop=False, tol=1e-12)[source]

Returns the hafnian of matrix with repeated rows/columns.

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:

>>> hafnian_repeated(A, rpt)
>>> hafnian(reduction(A, rpt))

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)

For more direct control, you may wish to call haf_rpt_real() or haf_rpt_complex() directly.

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.

  • use_eigen (bool) – if True (default), the Eigen linear algebra library is used for matrix multiplication. If the hafnian library was compiled with BLAS/Lapack support, then BLAS will be used for matrix multiplication.

  • tol (float) – the tolerance when checking that the matrix is symmetric. Default tolerance is 1e-12.

Returns

the hafnian of matrix A.

Return type

np.int64 or np.float64 or np.complex128

hafnian_batched(A, cutoff, mu=None, tol=1e-12, 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.

Returns

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

Return type

(array)

tor(A, fsum=False)[source]

Returns the Torontonian of a matrix.

For more direct control, you may wish to call tor_real() or tor_complex() directly.

The input matrix is cast to quadruple precision internally for a quadruple precision torontonian computation.

Parameters
  • A (array) – a np.complex128, square, symmetric array of even dimensions.

  • fsum (bool) – if True, the Shewchuck algorithm for more accurate summation is performed. This can significantly increase the accuracy of the computation, but no casting to quadruple precision takes place, as the Shewchuck algorithm only supports double precision.

Returns

the torontonian of matrix A.

Return type

np.float64 or np.complex128

perm(A, quad=True, fsum=False)[source]

Returns the permanent of a matrix via the Ryser formula.

For more direct control, you may wish to call perm_real() or perm_complex() directly.

Parameters
  • A (array) – a square array.

  • quad (bool) – If True, the input matrix is cast to a long double matrix internally for a quadruple precision hafnian computation.

  • fsum (bool) – Whether to use the fsum method for higher accuracy summation. Note that if fsum is true, double precision will be used, and the quad keyword argument will be ignored.

Returns

the permanent of matrix A.

Return type

np.float64 or np.complex128

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

np.int64 or np.float64 or np.complex128

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

hermite_multidimensional(R, cutoff, y=None, renorm=False, make_tensor=True, modified=False)[source]

Returns the multidimensional Hermite polynomials \(H_k^{(R)}(y)\).

Here \(R\) is an \(n \times n\) square matrix, and \(y\) is an \(n\) dimensional vector. The polynomials 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.

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

  • 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

Returns

the multidimensional Hermite polynomials

Return type

(array)

version()[source]

Get version number of The Walrus

Returns

The package version number

Return type

str