Source code for thewalrus._hafnian

# Copyright 2021 Xanadu Quantum Technologies Inc.

# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at

#     http://www.apache.org/licenses/LICENSE-2.0

# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
"""
Hafnian Python interface
"""
import warnings
from functools import lru_cache
from collections import Counter
from itertools import chain, combinations
import numba
import numpy as np
from thewalrus import charpoly


@numba.jit(nopython=True, cache=True)
def nb_binom(n, k):  # pragma: no cover
    """Numba version of binomial coefficient function.

    Args:
        n (int): how many options
        k (int): how many are chosen

    Returns:
        int: how many ways of choosing
    """
    if k < 0 or k > n:
        return 0
    if k in (0, n):
        return 1
    binom = 1
    for i in range(min(k, n - k)):
        binom *= n - i
        binom //= i + 1
    return binom


@numba.jit(nopython=True, cache=True)
def precompute_binoms(max_binom):  # pragma: no cover
    """Precompute binomial coefficients, return as a 2d array.

    Args:
        max_binom (int): max value of n in the binomial

    Returns:
        array: ``max_binom + 1 * max_binom + 1`` array of binomial coefficients
    """
    binoms = np.zeros((max_binom + 1, max_binom + 1), dtype=type(max_binom))
    for i in range(max_binom + 1):
        for j in range(max_binom + 1):
            binoms[i, j] = nb_binom(i, j)
    return binoms


@numba.jit(nopython=True, cache=True)
def nb_ix(arr, rows, cols):  # pragma: no cover
    """Numba implementation of ``np.ix_``.

    Args:
        arr (2d array): matrix to take submatrix of
        rows (array): rows to be selected in submatrix
        cols (array): columns to be selected in submatrix

    Return:
        array: selected submatrix of ``arr`` with dimension ``len(rows) * len(cols)``
    """
    return arr[rows][:, cols]


[docs]def matched_reps(reps): # pylint: disable = too-many-branches """Takes the repeated rows and find a way to pair them up to create a perfect matching with many repeated edges. Args: reps (list): list of repeated rows/cols Returns: tuple[array, array, int]: 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) """ n = len(reps) if sum(reps) == 0: return np.array([], dtype=int), np.array([], dtype=int), None # need to pair off the indices with high numbers of repetitions... x = range(n) # the starting set of indices edgesA = [] # contains part A of each pair edgesB = [] # part B of each pair edgereps = [] # number of repetitions of a pair reps, x = zip( *sorted(zip(reps, x), reverse=True) ) # sort according to reps, in descending order reps = list(reps) x = list(x) # remove zeros nonzero_reps = [] nonzero_x = [] for i, r in zip(x, reps): if r > 0: nonzero_reps.append(r) nonzero_x.append(i) reps = nonzero_reps x = nonzero_x while len(reps) > 1 or (len(reps) == 1 and reps[0] > 1): reps, x = zip(*sorted(zip(reps, x), reverse=True)) # sort reps = list(reps) x = list(x) if len(reps) == 1 or reps[0] > reps[1] * 2: # if largest number of reps is more than double the 2nd largest, pair it with itself edgesA += [x[0]] edgesB += [x[0]] edgereps += [reps[0] // 2] if reps[0] % 2 == 0: x = x[1:] reps = reps[1:] else: reps[0] = 1 else: # otherwise, form pairs between largest reps and 2nd largest reps edgesA += [x[0]] edgesB += [x[1]] edgereps += [reps[1]] if reps[0] > reps[1]: if len(x) > 2: x = [x[0]] + x[2:] reps = [reps[0] - reps[1]] + reps[2:] else: x = [x[0]] reps = [reps[0] - reps[1]] else: x = x[2:] reps = reps[2:] if len(x) == 1: oddmode = x[0] # if there is an unpaired mode, store it else: oddmode = None # the adjacency matrix of red edges connects 1 to N/2+1, 2 to N/2+2, etc. # Reorder the indices (from x2 back to x) so that the paired indices get # connected by red edges x = np.asarray(edgesA + edgesB, dtype=np.int64) # reordered list of indices edgereps = np.asarray(edgereps, dtype=np.int64) return x, edgereps, oddmode
@numba.jit(nopython=True, cache=True) def find_kept_edges(j, reps): # pragma: no cover """Write ``j`` as a string where the ith digit is in base ``reps[i]+1`` decides which edges are included given index of the inclusion/exclusion sum. Args: j (int): index of sum reps (list): number of repetitions of each edge Returns: array: number of repetitions kept for the current inclusion/exclusion step """ num = j output = [] bases = np.asarray(reps) + 1 for base in bases[::-1]: output.append(num % base) num //= base return np.array(output[::-1], dtype=reps.dtype) @numba.jit(nopython=True, cache=True) def f(A, n): # pragma: no cover """Evaluate the polynomial coefficients of the function in the eigenvalue-trace formula. Args: A (array): a two-dimensional matrix n (int): number of polynomial coefficients to compute Returns: array: polynomial coefficients """ # Compute combinations in O(n^2log n) time # code translated from thewalrus matlab script count = 0 comb = np.zeros((2, n // 2 + 1), dtype=np.complex128) comb[0, 0] = 1 powtrace = charpoly.powertrace(A, n // 2 + 1) for i in range(1, n // 2 + 1): factor = powtrace[i] / (2 * i) powfactor = 1 count = 1 - count comb[count, :] = comb[1 - count, :] for j in range(1, n // (2 * i) + 1): powfactor *= factor / j for k in range(i * j + 1, n // 2 + 2): comb[count, k - 1] += comb[1 - count, k - i * j - 1] * powfactor return comb[count, :] @numba.jit(nopython=True, cache=True) def f_loop(AX, AX_S, XD_S, D_S, n): # pragma: no cover """Evaluate the polynomial coefficients of the function in the eigenvalue-trace formula. Args: AX (array): two-dimensional matrix AX_S (array): ``AX_S`` with weights given by repetitions and excluded rows removed XD_S (array): diagonal multiplied by ``X`` D_S (array): diagonal n (int): number of polynomial coefficients to compute Returns: array: polynomial coefficients """ # Compute combinations in O(n^2log n) time # code translated from thewalrus matlab script count = 0 comb = np.zeros((2, n // 2 + 1), dtype=np.complex128) comb[0, 0] = 1 powtrace = charpoly.powertrace(AX, n // 2 + 1) for i in range(1, n // 2 + 1): factor = powtrace[i] / (2 * i) + (XD_S @ D_S) / 2 XD_S = XD_S @ AX_S powfactor = 1 count = 1 - count comb[count, :] = comb[1 - count, :] for j in range(1, n // (2 * i) + 1): powfactor *= factor / j for k in range(i * j + 1, n // 2 + 2): comb[count, k - 1] += comb[1 - count, k - i * j - 1] * powfactor return comb[count, :] # pylint: disable = too-many-arguments @numba.jit(nopython=True, cache=True) def f_loop_odd(AX, AX_S, XD_S, D_S, n, oddloop, oddVX_S): # pragma: no cover """Evaluate the polynomial coefficients of the function in the eigenvalue-trace formula when there is a self-edge in the fixed perfect matching. Args: AX (array): two-dimensional matrix AX_S (array): ``AX_S`` with weights given by repetitions and excluded rows removed XD_S (array): diagonal multiplied by ``X`` D_S (array): diagonal n (int): number of polynomial coefficients to compute oddloop (float): weight of self-edge oddVX_S (array): vector corresponding to matrix at the index of the self-edge Returns: array: polynomial coefficients """ count = 0 comb = np.zeros((2, n + 1), dtype=np.complex128) comb[0, 0] = 1 powtrace = charpoly.powertrace(AX, n + 1) for i in range(1, n + 1): if i == 1: factor = oddloop elif i % 2 == 0: factor = powtrace[i // 2] / i + (XD_S @ D_S) / 2 else: factor = oddVX_S @ D_S D_S = AX_S @ D_S powfactor = 1 count = 1 - count comb[count, :] = comb[1 - count, :] for j in range(1, n // i + 1): powfactor *= factor / j for k in range(i * j + 1, n + 2): comb[count, k - 1] += comb[1 - count, k - i * j - 1] * powfactor return comb[count, :] @numba.jit(nopython=True, cache=True) def get_AX_S(kept_edges, A): # pragma: no cover """Given the kept edges, return the appropriate scaled submatrices to compute ``f``. Args: kept_edges (array): number of repetitions of each edge A (array): matrix before repetitions applied Returns: array: scaled ``A @ X``, where ``X = ((0, I), (I, 0))`` """ z = np.concatenate((kept_edges, kept_edges)) nonzero_rows = np.where(z != 0)[0] n_nonzero_edges = len(nonzero_rows) // 2 kept_edges_nonzero = kept_edges[np.where(kept_edges != 0)] A_nonzero = nb_ix(A, nonzero_rows, nonzero_rows) AX_nonzero = np.empty_like(A_nonzero, dtype=np.complex128) AX_nonzero[:, :n_nonzero_edges] = kept_edges_nonzero * A_nonzero[:, n_nonzero_edges:] AX_nonzero[:, n_nonzero_edges:] = kept_edges_nonzero * A_nonzero[:, :n_nonzero_edges] return AX_nonzero @numba.jit(nopython=True, cache=True) def get_submatrices(kept_edges, A, D, oddV): # pragma: no cover """Given the kept edges, return the appropriate scaled submatrices to compute ``f``. Args: kept_edges (array): number of repetitions of each edge A (array): matrix before repetitions applied D (array): diagonal before repetitions applied oddV (array): Row of matrix at index of self-edge. ``None`` is no self-edge. Returns: tuple[array, array, array, array]: scaled ``A @ X `` (where ``X = ((0, I), (I, 0))``), scaled ``X @ D``, scaled ``D``, and scaled ``oddV @ X`` """ z = np.concatenate((kept_edges, kept_edges)) nonzero_rows = np.where(z != 0)[0] n_nonzero_edges = len(nonzero_rows) // 2 kept_edges_nonzero = kept_edges[np.where(kept_edges != 0)] A_nonzero = nb_ix(A, nonzero_rows, nonzero_rows) AX_nonzero = np.empty_like(A_nonzero, dtype=np.complex128) AX_nonzero[:, :n_nonzero_edges] = kept_edges_nonzero * A_nonzero[:, n_nonzero_edges:] AX_nonzero[:, n_nonzero_edges:] = kept_edges_nonzero * A_nonzero[:, :n_nonzero_edges] D_nonzero = D[nonzero_rows] XD_nonzero = np.empty_like(D_nonzero, dtype=np.complex128) XD_nonzero[:n_nonzero_edges] = kept_edges_nonzero * D_nonzero[n_nonzero_edges:] XD_nonzero[n_nonzero_edges:] = kept_edges_nonzero * D_nonzero[:n_nonzero_edges] if oddV is not None: oddV_nonzero = oddV[nonzero_rows] oddVX_nonzero = np.empty_like(oddV_nonzero, dtype=np.complex128) oddVX_nonzero[:n_nonzero_edges] = kept_edges_nonzero * oddV_nonzero[n_nonzero_edges:] oddVX_nonzero[n_nonzero_edges:] = kept_edges_nonzero * oddV_nonzero[:n_nonzero_edges] else: oddVX_nonzero = None return AX_nonzero, XD_nonzero, D_nonzero, oddVX_nonzero @numba.jit(nopython=True, cache=True) def get_submatrix_batch_odd0(kept_edges, oddV0): # pragma: no cover """Find ``oddVX_nonzero0`` for batching (sometimes different vertices are identified as self edges). Args: kept_edges (array): number of repetitions of each edge oddV0 (array): Row of matrix at index of self-edge. ``None`` is no self-edge. Returns: array: scaled ``oddV0 @ X`` """ z = np.concatenate((kept_edges, kept_edges)) nonzero_rows = np.where(z != 0)[0] n_nonzero_edges = len(nonzero_rows) // 2 kept_edges_nonzero = kept_edges[np.where(kept_edges != 0)] oddV_nonzero0 = oddV0[nonzero_rows] oddVX_nonzero0 = np.empty_like(oddV_nonzero0, dtype=np.complex128) oddVX_nonzero0[:n_nonzero_edges] = kept_edges_nonzero * oddV_nonzero0[n_nonzero_edges:] oddVX_nonzero0[n_nonzero_edges:] = kept_edges_nonzero * oddV_nonzero0[:n_nonzero_edges] return oddVX_nonzero0 @numba.jit(nopython=True, cache=True) def get_Dsubmatrices(kept_edges, D): # pragma: no cover """Find submatrices for batch gamma functions.""" z = np.concatenate((kept_edges, kept_edges)) nonzero_rows = np.where(z != 0)[0] n_nonzero_edges = len(nonzero_rows) // 2 kept_edges_nonzero = kept_edges[np.where(kept_edges != 0)] D_nonzero = D[nonzero_rows] XD_nonzero = np.empty_like(D_nonzero, dtype=np.complex128) XD_nonzero[:n_nonzero_edges] = kept_edges_nonzero * D_nonzero[n_nonzero_edges:] XD_nonzero[n_nonzero_edges:] = kept_edges_nonzero * D_nonzero[:n_nonzero_edges] return XD_nonzero, D_nonzero @numba.jit(nopython=True, cache=True) def eigvals(M): # pragma: no cover """Computes the eigenvalues of a matrix. Args: M (array): square matrix Returns: array: eigenvalues of the matrix ``M`` """ return np.linalg.eigvals(M) # pylint: disable=W0612, E1133 @numba.jit(nopython=True, parallel=True, cache=True) def _calc_hafnian(A, edge_reps, glynn=True): # pragma: no cover r"""Compute hafnian, using inputs as prepared by frontend hafnian function compiled with Numba. Args: A (array): matrix ordered according to the chosen perfect matching edge_reps (array): how many times each edge in the perfect matching is repeated glynn (bool): whether to use finite difference sieve Returns: complex: value of hafnian """ n = A.shape[0] N = 2 * edge_reps.sum() # number of photons if glynn: steps = ((edge_reps[0] + 2) // 2) * np.prod(edge_reps[1:] + 1) else: steps = np.prod(edge_reps + 1) # precompute binomial coefficients max_binom = edge_reps.max() + 1 binoms = precompute_binoms(max_binom) H = np.complex128(0) # start running total for the hafnian for j in numba.prange(steps): kept_edges = find_kept_edges(j, edge_reps) edge_sum = kept_edges.sum() binom_prod = 1.0 for i in range(n // 2): binom_prod *= binoms[edge_reps[i], kept_edges[i]] if glynn: kept_edges = 2 * kept_edges - edge_reps AX_S = get_AX_S(kept_edges, A) prefac = (-1.0) ** (N // 2 - edge_sum) * binom_prod if glynn and kept_edges[0] == 0: prefac *= 0.5 Hnew = prefac * f(AX_S, N)[N // 2] H += Hnew if glynn: H = H * 0.5 ** (N // 2 - 1) return H def _haf(A, reps=None, glynn=True): r"""Calculate hafnian with (optional) repeated rows and columns. Code contributed by `Jake F.F. Bulmer <https://github.com/jakeffbulmer/gbs>`_ based on `arXiv:2108.01622 <https://arxiv.org/abs/2108.01622>`_. Args: A (array): N x N matrix. 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 hafnian calculation """ n = A.shape[0] if reps is None: reps = [1] * n N = sum(reps) if N == 0: return 1.0 if N % 2 == 1: return 0.0 assert n == len(reps) x, edge_reps, oddmode = matched_reps(reps) # make new A matrix using the ordering from above Ax = A[np.ix_(x, x)].astype(np.complex128) H = _calc_hafnian(Ax, edge_reps, glynn) return H # pylint: disable=too-many-arguments, redefined-outer-name, not-an-iterable @numba.jit(nopython=True, parallel=True, cache=True) def _calc_loop_hafnian(A, D, edge_reps, oddloop=None, oddV=None, glynn=True): # pragma: no cover """Compute loop hafnian, using inputs as prepared by frontend loop_hafnian function compiled with Numba. Code contributed by `Jake F.F. Bulmer <https://github.com/jakeffbulmer/gbs>`_ based on `arXiv:2108.01622 <https://arxiv.org/abs/2108.01622>`_. Args: A (array): matrix ordered according to the chosen perfect matching. D (array): diagonals ordered according to the chosen perfect matchin edge_reps (array): how many times each edge in the perfect matching is repeated oddloop (float): weight of self-loop in perfect matching, None if no self-loops oddV (array): row of matrix corresponding to the odd loop in the perfect matching glynn (bool): whether to use finite difference sieve Returns: complex: value of loop hafnian """ n = A.shape[0] N = 2 * edge_reps.sum() # Number of photons if oddloop is not None: N += 1 if glynn and (oddloop is None): steps = ((edge_reps[0] + 2) // 2) * np.prod(edge_reps[1:] + 1) else: steps = np.prod(edge_reps + 1) # Precompute binomial coefficients max_binom = edge_reps.max() + 1 binoms = precompute_binoms(max_binom) H = np.complex128(0) # Start running total for the hafnian for j in numba.prange(steps): kept_edges = find_kept_edges(j, edge_reps) edge_sum = kept_edges.sum() binom_prod = 1.0 for i in range(n // 2): binom_prod *= binoms[edge_reps[i], kept_edges[i]] if glynn: kept_edges = 2 * kept_edges - edge_reps AX_S, XD_S, D_S, oddVX_S = get_submatrices(kept_edges, A, D, oddV) AX = AX_S.copy() prefac = (-1.0) ** (N // 2 - edge_sum) * binom_prod if oddloop is not None: Hnew = prefac * f_loop_odd(AX, AX_S, XD_S, D_S, N, oddloop, oddVX_S)[N] else: if glynn and kept_edges[0] == 0: prefac *= 0.5 Hnew = prefac * f_loop(AX, AX_S, XD_S, D_S, N)[N // 2] H += Hnew if glynn: if oddloop is None: H = H * 0.5 ** (N // 2 - 1) else: H = H * 0.5 ** (N // 2) return H # pylint: disable=redefined-outer-name
[docs]def loop_hafnian(A, D=None, reps=None, glynn=True): """Calculate loop hafnian with (optional) repeated rows and columns. Code contributed by `Jake F.F. Bulmer <https://github.com/jakeffbulmer/gbs>`_ based on `arXiv:2108.01622 <https://arxiv.org/abs/2108.01622>`_. Args: 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 """ n = A.shape[0] if reps is None: reps = [1] * n if D is None: D = A.diagonal() N = sum(reps) if N == 0: return 1.0 if N == 1: return D[np.where(np.array(reps) == 1)[0][0]] assert n == len(reps) assert D.shape[0] == n x, edge_reps, oddmode = matched_reps(reps) # Make new A matrix and D vector using the ordering from above if oddmode is not None: oddloop = D[oddmode].astype(np.complex128) oddV = A[oddmode, x].astype(np.complex128) else: oddloop = None oddV = None Ax = A[np.ix_(x, x)].astype(np.complex128) Dx = D[x].astype(np.complex128) H = _calc_loop_hafnian(Ax, Dx, edge_reps, oddloop, oddV, glynn) return H
def input_validation(A, rtol=1e-05, atol=1e-08): """Checks that the matrix A satisfies the requirements for Hafnian calculation. These include: * That the ``A`` is a NumPy array * That ``A`` is square * That ``A`` does not contain any NaNs * That ``A`` is symmetric Args: A (array): a NumPy array. rtol (float): the relative tolerance parameter used in ``np.allclose`` atol (float): the absolute tolerance parameter used in ``np.allclose`` Returns: bool: returns ``True`` if the matrix satisfies all requirements """ if not isinstance(A, np.ndarray): raise TypeError("Input matrix must be a NumPy array.") n = A.shape if n[0] != n[1]: raise ValueError("Input matrix must be square.") if np.isnan(A).any(): raise ValueError("Input matrix must not contain NaNs.") if not np.allclose(A, A.T, rtol=rtol, atol=atol): raise ValueError("Input matrix must be symmetric.") return True def bandwidth(A): """Calculates the upper bandwidth of the matrix A. Args: A (array): input matrix Returns: int: bandwidth of matrix """ n, _ = A.shape for i in range(n - 1, 0, -1): vali = np.diag(A, i) if not np.allclose(vali, 0): return i return 0 def powerset(iterable): """Calculates the powerset of a list. Args: iterable (iterable): input list Returns: chain: chain of all subsets of input list """ return chain.from_iterable(combinations(iterable, r) for r in range(len(iterable) + 1))
[docs]def reduction(A, rpt): r"""Calculates the reduction of an array by a vector of indices. This is equivalent to repeating the ith row/column of :math:`A`, :math:`rpt_i` times. Args: 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: array: the reduction of ``A`` by the index vector ``rpt`` """ rows = [i for sublist in [[idx] * j for idx, j in enumerate(rpt)] for i in sublist] if A.ndim == 1: return A[rows] return A[:, rows][rows]
# pylint: disable=too-many-arguments
[docs]def hafnian( A, loop=False, rtol=1e-05, atol=1e-08, approx=False, num_samples=1000, method="glynn", ): # pylint: disable=too-many-arguments """Returns the hafnian of a matrix. Code contributed by `Jake F.F. Bulmer <https://github.com/jakeffbulmer/gbs>`_ based on `arXiv:2108.01622 <https://arxiv.org/abs/2108.01622>`_. Args: 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: int or float or complex: the hafnian of matrix ``A`` """ # pylint: disable=too-many-return-statements,too-many-branches input_validation(A, rtol=rtol, atol=atol) matshape = A.shape if method == "glynn": glynn = True if method == "inclexcl": glynn = False if matshape == (0, 0): return 1 if matshape[0] % 2 != 0 and not loop: return 0.0 if np.allclose(np.diag(np.diag(A)), A, rtol=rtol, atol=atol): if loop: return np.prod(np.diag(A)) return 0 matshape = A.shape if matshape[0] == 2: if loop: return A[0, 1] + A[0, 0] * A[1, 1] return A[0][1] if matshape[0] == 3 and loop: return ( A[0, 0] * A[1, 2] + A[1, 1] * A[0, 2] + A[2, 2] * A[0, 1] + A[0, 0] * A[1, 1] * A[2, 2] ) if matshape[0] == 4: if loop: result = ( A[0, 1] * A[2, 3] + A[0, 2] * A[1, 3] + A[0, 3] * A[1, 2] + A[0, 0] * A[1, 1] * A[2, 3] + A[0, 1] * A[2, 2] * A[3, 3] + A[0, 2] * A[1, 1] * A[3, 3] + A[0, 0] * A[2, 2] * A[1, 3] + A[0, 0] * A[3, 3] * A[1, 2] + A[0, 3] * A[1, 1] * A[2, 2] + A[0, 0] * A[1, 1] * A[2, 2] * A[3, 3] ) return result return A[0, 1] * A[2, 3] + A[0, 2] * A[1, 3] + A[0, 3] * A[1, 2] if approx: if np.any(np.iscomplex(A)): raise ValueError("Input matrix must be real") if np.any(A < 0): raise ValueError("Input matrix must not have negative entries") return hafnian_approx(A, num_samples=num_samples) if loop: if method == "recursive": warnings.warn("Recursive algorithm does not support the loop hafnian") return loop_hafnian(A, D=None, reps=None, glynn=True) if method == "recursive": return recursive_hafnian(A) return _haf(A, reps=None, glynn=glynn)
[docs]def hafnian_sparse(A, D=None, loop=False): r"""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 :func:`~.hafnian` happens around 50% sparsity. Args: 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: float: hafnian of ``A`` or of the submatrix of ``A`` defined by the set of indices ``D`` """ if D is None: D = frozenset(range(len(A))) else: D = frozenset(D) if not loop: A = A - np.diag(np.diag(A)) if np.allclose(A, 0): return 0.0 r, _ = np.nonzero(A) m = max(Counter(r).values()) # max nonzero values per row/column @lru_cache(maxsize=2**m) def indices(d, k): return d.intersection(set(np.nonzero(A[k])[0])) @lru_cache(maxsize=2**m) def lhaf(d: frozenset) -> float: if not d: return 1 d_ = set(d) k = d_.pop() return sum(A[i, k] * lhaf(frozenset(d_).difference({i})) for i in indices(d, k)) return lhaf(D)
[docs]def hafnian_repeated(A, rpt, mu=None, loop=False, rtol=1e-05, atol=1e-08, glynn=True): r"""Returns the hafnian of matrix with repeated rows/columns. Code contributed by `Jake F.F. Bulmer <https://github.com/jakeffbulmer/gbs>`_ based on `arXiv:2108.01622 <https://arxiv.org/abs/2108.01622>`_. The :func:`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: .. code: >>> 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 (:math:`\sum_{i}rpt_i \gg N`) can be significantly faster. .. note:: If :math:`rpt=(1, 1, \dots, 1)`, then >>> hafnian_repeated(A, rpt) == hafnian(A) Args: A (array): a square, symmetric :math:`N\times N` array rpt (Sequence): a length-:math:`N` positive integer sequence, corresponding to the number of times each row/column of matrix :math:`A` is repeated mu (array): A vector of length :math:`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: int or float or complex: the hafnian of matrix A """ # pylint: disable=too-many-return-statements,too-many-branches input_validation(A, atol=atol, rtol=rtol) if len(rpt) != len(A): raise ValueError("the rpt argument must be 1-dimensional sequence of length len(A).") nud = np.array(rpt, dtype=np.int32) if not np.all(np.mod(rpt, 1) == 0) or np.any(nud < 0): raise ValueError("the rpt argument must contain non-negative integers.") if np.all(nud == 0): return 1.0 if np.sum(nud) % 2 != 0 and not loop: return 0.0 if mu is None: mu = A.diagonal().copy() if np.allclose(A, 0, rtol=rtol, atol=atol): if loop: return np.prod(mu**rpt) return 0 if len(mu) != len(A): raise ValueError("Length of means vector must be the same length as the matrix A.") if loop: return loop_hafnian(A, D=mu, reps=rpt, glynn=glynn) return _haf(A, reps=rpt, glynn=glynn)
[docs]def hafnian_banded(A, loop=False, rtol=1e-05, atol=1e-08): """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. <https://arxiv.org/abs/2009.11824>`_. Args: A (array): a square, symmetric array of even dimensions Returns: int or float or complex: the loop hafnian of matrix ``A`` """ input_validation(A, atol=atol, rtol=rtol) (n, _) = A.shape w = bandwidth(A) if not loop: A = A - np.diag(np.diag(A)) loop_haf = {(): 1, (1,): A[0, 0]} for t in range(1, n + 1): if t - 2 * w - 1 > 0: lower_end = set(range(1, t - 2 * w)) else: lower_end = set() upper_end = set(range(1, t + 1)) diff = [item for item in upper_end if item not in lower_end] # Makes sure set ordering is preserved when the difference of two sets is taken # This is also used in the if statement below ps = powerset(diff) lower_end = tuple(lower_end) for D in ps: if lower_end + D not in loop_haf: loop_haf[lower_end + D] = sum( [ A[i - 1, t - 1] * loop_haf[tuple(item for item in lower_end + D if item not in set((i, t)))] for i in D ] ) return loop_haf[tuple(range(1, n + 1))]
@numba.jit(nopython=True) def recursive_hafnian(A): # pragma: no cover r"""Computes the hafnian of the matrix with the recursive algorithm. It is an implementation of algorithm 2 in *Counting perfect matchings as fast as Ryser* :cite:`bjorklund2012counting`. This code is a modified version of the code found here: `Recursive hafnian <https://codegolf.stackexchange.com/questions/157049/calculate-the-hafnian-as-quickly-as-possible>`_. Args: A (array): the input matrix Returns: float: the hafnian of the input matrix """ nb_lines, nb_columns = A.shape if nb_lines != nb_columns: raise ValueError("Matrix must be square") if nb_lines % 2 != 0: raise ValueError("Matrix size must be even") n = len(A) // 2 z = np.zeros((n * (2 * n - 1), n + 1), dtype=A.dtype) for j in range(1, 2 * n): ind = j * (j - 1) // 2 for k in range(j): z[ind + k][0] = A[j][k] g = np.zeros(n + 1, dtype=A.dtype) g[0] = 1 return solve(z, 2 * n, 1, g, n) @numba.jit(nopython=True) def solve(b, s, w, g, n): # pragma: no cover r"""Implements the recursive algorithm. Args: b (array): matrix that is transformed recursively s (int): size of the original matrix that changes at every recursion k (int): a variable of the recursive algorithm g (int): matrix that is transformed recursively n (int): size of the original matrix divided by 2 Returns: float: the hafnian of the input matrix """ if s == 0: return w * g[n] c = np.zeros(((s - 2) * (s - 3) // 2, n + 1), dtype=g.dtype) i = 0 for j in range(1, s - 2): for k in range(j): c[i] = b[(j + 1) * (j + 2) // 2 + k + 2] i += 1 h = solve(c, s - 2, -w, g, n) e = g.copy() for u in range(n): for v in range(n - u): e[u + v + 1] += g[u] * b[0][v] for j in range(1, s - 2): for k in range(j): for u in range(n): for v in range(n - u): c[j * (j - 1) // 2 + k][u + v + 1] += ( b[(j + 1) * (j + 2) // 2][u] * b[(k + 1) * (k + 2) // 2 + 1][v] + b[(k + 1) * (k + 2) // 2][u] * b[(j + 1) * (j + 2) // 2 + 1][v] ) return h + solve(c, s - 2, w, e, n) @numba.jit(nopython=True) def _one_det(B): # pragma: no cover """Calculates the determinant of an antisymmetric matrix with entries distributed according to a normal distribution, with scale equal to the entries of the symmetric matrix given as input. Args: B (array[float]): symmetric matrix Returns: float: determinant of the samples antisymmetric matrix """ mat = np.empty_like(B, dtype=np.float64) n, m = B.shape for i in range(n): for j in range(m): mat[i, j] = B[i, j] * np.random.normal() mat[j, i] = -mat[i, j] return np.linalg.det(mat) @numba.jit(nopython=True) def hafnian_approx(A, num_samples=1000): # pragma: no cover """Returns the approximation to the hafnian of a matrix with non-negative entries. The approximation follows the stochastic Barvinok's approximation allowing the hafnian can be approximated as the sum of determinants of matrices. The accuracy of the approximation increases with increasing number of iterations. Args: B (array[float]): a symmetric matrix Returns: float: approximate hafnian of the input """ sqrtA = np.sqrt(A) return np.array([_one_det(sqrtA) for _ in range(num_samples)]).mean()