search

XanaduAI / thewalrus

A library for the calculation of hafnians, Hermite polynomials and Gaussian boson sampling.
https://the-walrus.readthedocs.io
Apache License 2.0
99 stars 54 forks source link

RuntimeWarning: divide by zero encountered in det #367

Closed minseok1999 closed 1 year ago

minseok1999 commented 1 year ago

Before posting a bug report

Expected behavior

I first prepared arbitrary adjacency matrix of a graph as follows.

matrix = np.array( [[1,1,1,1,1,1,0,0,0,0], [1,1,1,1,1,1,0,0,0,0], [1,1,1,1,1,1,0,0,0,0], [1,1,1,1,1,1,0,0,1,0], [1,1,1,1,1,1,0,1,0,0], [1,1,1,1,1,1,1,0,0,0], [0,0,0,0,0,1,1,0,0,0], [0,0,0,0,1,0,0,1,0,0], [0,0,0,1,0,0,0,0,1,1], [0,0,0,0,0,0,0,0,1,1]] )-np.eye(10)

And then I fed this matrix into the module

s=hafnian_sample_graph(matrix,6) print(s)

to see the simulated result of photo count with mean photon number of 6 at the output Gaussian state.

I have only changed the number of sample parameter in the module as follows to generate 10 samples

def hafnian_sample_graph( A, n_mean, samples=10, cutoff=5, max_photons=30, approx=False, approx_samples=1e5, pool=False ):

Actual behavior

But Implementing this line gave this error message

/Users/ryuminseok/anaconda3/lib/python3.10/site-packages/numpy/linalg/linalg.py:2154: RuntimeWarning: divide by zero encountered in det r = _umath_linalg.det(a, signature=signature) /Users/ryuminseok/anaconda3/lib/python3.10/site-packages/numpy/linalg/linalg.py:2154: RuntimeWarning: invalid value encountered in det r = _umath_linalg.det(a, signature=signature) /Users/ryuminseok/anaconda3/lib/python3.10/site-packages/numpy/linalg/linalg.py:2154: RuntimeWarning: divide by zero encountered in det r = _umath_linalg.det(a, signature=signature) /Users/ryuminseok/anaconda3/lib/python3.10/site-packages/numpy/linalg/linalg.py:2154: RuntimeWarning: invalid value encountered in det r = _umath_linalg.det(a, signature=signature)

The code still generates output pattern, although I have to wait for quite a long time. But I doubt that this module of photo count pattern generation is reliable if it involves division by zero which can lead to numerical instability or incorrect result.

Reproduces how often

It occurs every time

System information

The Walrus: a Python library for for the calculation of hafnians, Hermite polynomials, and Gaussian boson sampling. Copyright 2018-2021 Xanadu Quantum Technologies Inc.

Python version: 3.10.9 Platform info: macOS-13.4.1-arm64-arm-64bit Installation path: /Users/ryuminseok/anaconda3/lib/python3.10/site-packages/thewalrus The Walrus version: 0.20.0 Numpy version: 1.23.5 Scipy version: 1.10.0 SymPy version: 1.11.1 Numba version: 0.56.4

Source code

No response

Tracebacks

No response

Additional information

I have also tried out first generating covariance matrix corresponding to the adjacency matrix of my interest by 

Q=gen_Qmat_from_graph(matrix,6)
cov=Covmat(Q)
sample=hafnian_sample_state(cov,1)
print(sample)

But it gave same error message
isaacdevlugt commented 1 year ago

Hey @minseok1999! Can you please edit your post and include the output of the following:

import thewalrus as tw
tw.about()

You can put this in the "System information" box. It would also help if you could include the full error traceback 😄.

minseok1999 commented 1 year ago

Hey @minseok1999! Can you please edit your post and include the output of the following:

import thewalrus as tw
tw.about()

You can put this in the "System information" box. It would also help if you could include the full error traceback 😄.

I have edited as requested!

isaacdevlugt commented 1 year ago

Hmm... I'm not able to replicate the behaviour you're seeing. Here are my package versions:

Python version:            3.9.0
The Walrus version:        0.20.0
Numpy version:             1.23.5
Scipy version:             1.10.1
SymPy version:             1.12
Numba version:             0.57.0

Slightly different from yours. Can you try those versions and see if it works for you?

minseok1999 commented 1 year ago

Hmm... I'm not able to replicate the behaviour you're seeing. Here are my package versions:

Python version:            3.9.0
The Walrus version:        0.20.0
Numpy version:             1.23.5
Scipy version:             1.10.1
SymPy version:             1.12
Numba version:             0.57.0

Slightly different from yours. Can you try those versions and see if it works for you?

Here is the list of libraries that I imported

import multiprocessing
from multiprocessing import Pool

import numpy as np
from scipy.special import factorial as fac

from thewalrus._hafnian import hafnian, reduction
from thewalrus._torontonian import tor
from thewalrus.quantum import (
    Amat,
    Covmat,
    Qmat,
    Xmat,
    gen_Qmat_from_graph,
    is_classical_cov,
    reduced_gaussian,
    density_matrix_element,
)

# ===============================================================================================
# Hafnian sampling
# ===============================================================================================

# pylint: disable=too-many-branches
def generate_hafnian_sample(
    cov, mean, hbar=2, cutoff=6, max_photons=30, approx=False, approx_samples=1e5
):  # pylint: disable=too-many-branches
    r"""Returns a single sample from the Hafnian of a Gaussian state.

    Args:
        cov (array): a :math:`2N\times 2N` ``np.float64`` covariance matrix
            representing an :math:`N` mode quantum state. This can be obtained
            via the ``scovmavxp`` method of the Gaussian backend of Strawberry Fields.
        mean (array): a :math:`2N`` ``np.float64`` vector of means representing the Gaussian
            state.
        hbar (float): (default 2) the value of :math:`\hbar` in the commutation
            relation :math:`[\x,\p]=i\hbar`.
        cutoff (int): the Fock basis truncation.
        max_photons (int): specifies the maximum number of photons that can be counted.
        approx (bool): if ``True``, the approximate hafnian algorithm is used.
            Note that this can only be used for real, non-negative matrices.
        approx_samples: the number of samples used to approximate the hafnian if ``approx=True``.

    Returns:
        np.array[int]: a photon number sample from the Gaussian states.
    """
    N = len(cov) // 2
    result = []
    prev_prob = 1.0
    nmodes = N
    if mean is None:
        local_mu = np.zeros(2 * N)
    else:
        local_mu = mean
    A = Amat(Qmat(cov), hbar=hbar)

    for k in range(nmodes):
        probs1 = np.zeros([cutoff + 1], dtype=np.float64)
        kk = np.arange(k + 1)
        mu_red, V_red = reduced_gaussian(local_mu, cov, kk)

        if approx:
            Q = Qmat(V_red, hbar=hbar)
            A = Amat(Q, hbar=hbar, cov_is_qmat=True)

        for i in range(cutoff):
            indices = result + [i]
            ind2 = indices + indices
            if approx:
                factpref = np.prod(fac(indices))
                mat = reduction(A, ind2)
                probs1[i] = (
                    hafnian(np.abs(mat.real), approx=True, num_samples=approx_samples) / factpref
                )
            else:
                probs1[i] = density_matrix_element(
                    mu_red, V_red, indices, indices, include_prefactor=True, hbar=hbar
                ).real

        if approx:
            probs1 = probs1 / np.sqrt(np.linalg.det(Q).real)

        probs2 = probs1 / prev_prob
        probs3 = np.maximum(
            probs2, np.zeros_like(probs2)
        )  # pylint: disable=assignment-from-no-return
        ssum = np.sum(probs3)
        if ssum < 1.0:
            probs3[-1] = 1.0 - ssum

        # The following normalization of probabilities is needed when approx=True
        if approx:
            if ssum > 1.0:
                probs3 = probs3 / ssum

        result.append(np.random.choice(a=range(len(probs3)), p=probs3))
        if result[-1] == cutoff:
            return -1
        if sum(result) > max_photons:
            return -1
        prev_prob = probs1[result[-1]]

    return result

def _hafnian_sample(args):
    r"""Returns samples from the Hafnian of a Gaussian state.

    Note: this is a wrapper function, instead of using this function
    directly, please use either :func:`torontonian_sample_state` or
    :func:`torontonian_sample_graph`.

    Args:
        args (list): a list containing the following parameters:

            cov (array)
                a :math:`2N\times 2N` ``np.float64`` covariance matrix
                representing an :math:`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 :math:`2N`` ``np.float64`` vector of means representing the Gaussian
                state.

            hbar (float)
                the value of :math:`\hbar` in the commutation relation :math:`[\x,\p]=i\hbar`.

            cutoff (int)
                the Fock basis truncation.

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

            approx (bool)
                if ``True``, the approximate hafnian algorithm is used.
                Note that this can only be used for real, non-negative matrices.

            approx_samples (int)
                the number of samples used to approximate the hafnian if ``approx=True``.

    Returns:
        np.array[int]: photon number samples from the Gaussian state
    """
    cov, samples, mean, hbar, cutoff, max_photons, approx, approx_samples = args

    if not isinstance(cov, np.ndarray):
        raise TypeError("Covariance matrix must be a NumPy array.")

    matshape = cov.shape

    if matshape[0] != matshape[1]:
        raise ValueError("Covariance matrix must be square.")

    if np.isnan(cov).any():
        raise ValueError("Covariance matrix must not contain NaNs.")

    samples_array = []
    j = 0

    while j < samples:
        result = generate_hafnian_sample(
            cov,
            mean=mean,
            hbar=hbar,
            cutoff=cutoff,
            max_photons=max_photons,
            approx=approx,
            approx_samples=approx_samples,
        )

        if result != -1:
            # if result == -1, then you never get anything beyond cutoff
            samples_array.append(result)
            j = j + 1

    return np.vstack(samples_array)

def hafnian_sample_state(
    cov,
    samples,
    mean=None,
    hbar=2,
    cutoff=5,
    max_photons=30,
    approx=False,
    approx_samples=1e5,
    pool=False,
):
    r"""Returns samples from the Hafnian of a Gaussian state.

    Args:
        cov (array): a :math:`2N\times 2N` ``np.float64`` covariance matrix
            representing an :math:`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 :math:`2N`` ``np.float64`` vector of means representing the Gaussian
                state.
        hbar (float): (default 2) the value of :math:`\hbar` in the commutation
            relation :math:`[\x,\p]=i\hbar`.
        cutoff (int): the Fock basis truncation.
        max_photons (int): specifies the maximum number of photons that can be counted.
        approx (bool): if ``True``, the :func:`~.hafnian_approx` function is used
            to approximate the hafnian. Note that this can only be used for
            real, non-negative matrices.
        approx_samples: the number of samples used to approximate the hafnian if ``approx=True``.
        pool (bool): if ``True``, uses ``multiprocessor.Pool`` for parallelization of samples

    Returns:
        np.array[int]: photon number samples from the Gaussian state
    """
    if not pool:
        params = [cov, samples, mean, hbar, cutoff, max_photons, approx, approx_samples]
        return _hafnian_sample(params)

    pool = Pool()
    nprocs = multiprocessing.cpu_count()
    localsamps = samples // nprocs

    params = [[cov, localsamps, mean, hbar, cutoff, max_photons, approx, approx_samples]] * (nprocs - 1)
    params.append(
        [
            cov,
            samples - localsamps * (nprocs - 1),
            mean,
            hbar,
            cutoff,
            max_photons,
            approx,
            approx_samples,
        ]
    )

    result = np.vstack(pool.map(_hafnian_sample, params))
    pool.close()  # no more tasks
    pool.join()  # wrap up current tasks

    return result

def hafnian_sample_graph(
    A, n_mean, samples=10, cutoff=5, max_photons=30, approx=False, approx_samples=1e5, pool=False
):
    r"""Returns samples from the Gaussian state specified by the adjacency matrix :math:`A`
    and with total mean photon number :math:`n_{mean}`

    Args:
        A (array): a :math:`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.
        approx (bool): if ``True``, the approximate hafnian algorithm is used.
            Note that this can only be used for real, non-negative matrices.
        approx_samples: the number of samples used to approximate the hafnian if ``approx=True``.
        pool (bool): if ``True``, uses ``multiprocessor.Pool`` for parallelization of samples

    Returns:
        np.array[int]: photon number samples from the Gaussian state
    """
    Q = gen_Qmat_from_graph(A, n_mean)
    cov = Covmat(Q)
    return hafnian_sample_state(
        cov,
        samples,
        mean=None,
        hbar=2,
        cutoff=cutoff,
        max_photons=max_photons,
        approx=approx,
        approx_samples=approx_samples,
        pool=pool,
    )

def seed(seed_val=None):
    r""" Seeds the random number generator used in the sampling algorithms.

    This function is a wrapper around ``numpy.random.seed()``. By setting the seed
    to a specific integer, the sampling algorithms will exhibit deterministic behaviour.

    Args:
        seed_val (int): Seed for RandomState. Must be convertible to 32 bit unsigned integers.
    """
    np.random.seed(seed_val)

and I have changed the package version as follows

The Walrus: a Python library for for the calculation of hafnians, Hermite polynomials, and Gaussian boson sampling. Copyright 2018-2021 Xanadu Quantum Technologies Inc.

Python version: 3.10.9 Platform info: macOS-13.4.1-arm64-arm-64bit Installation path: /Users/ryuminseok/anaconda3/lib/python3.10/site-packages/thewalrus The Walrus version: 0.20.0 Numpy version: 1.23.5 Scipy version: 1.11.0 SymPy version: 1.12 Numba version: 0.57.1

But I keep seeing this behavior

/Users/ryuminseok/anaconda3/lib/python3.10/site-packages/numpy/linalg/linalg.py:2154: RuntimeWarning: divide by zero encountered in det
  r = _umath_linalg.det(a, signature=signature)
/Users/ryuminseok/anaconda3/lib/python3.10/site-packages/numpy/linalg/linalg.py:2154: RuntimeWarning: invalid value encountered in det
  r = _umath_linalg.det(a, signature=signature)
OMP: Info #276: omp_set_nested routine deprecated, please use omp_set_max_active_levels instead.
/Users/ryuminseok/anaconda3/lib/python3.10/site-packages/numpy/linalg/linalg.py:2154: RuntimeWarning: divide by zero encountered in det
  r = _umath_linalg.det(a, signature=signature)
/Users/ryuminseok/anaconda3/lib/python3.10/site-packages/numpy/linalg/linalg.py:2154: RuntimeWarning: invalid value encountered in det
  r = _umath_linalg.det(a, signature=signature)

Can this be because of the specific property of the matrix that I fed into hafnium_sample_graph as s=hafnian_sample_graph(matrix,6) ?

isaacdevlugt commented 1 year ago

@minseok1999 Let's stick with the simpler example here:

import numpy as np
from thewalrus.samples import hafnian_sample_graph

matrix = np.array(
    [
        [1, 1, 1, 1, 1, 1, 0, 0, 0, 0],
        [1, 1, 1, 1, 1, 1, 0, 0, 0, 0],
        [1, 1, 1, 1, 1, 1, 0, 0, 0, 0],
        [1, 1, 1, 1, 1, 1, 0, 0, 1, 0],
        [1, 1, 1, 1, 1, 1, 0, 1, 0, 0],
        [1, 1, 1, 1, 1, 1, 1, 0, 0, 0],
        [0, 0, 0, 0, 0, 1, 1, 0, 0, 0],
        [0, 0, 0, 0, 1, 0, 0, 1, 0, 0],
        [0, 0, 0, 1, 0, 0, 0, 0, 1, 1],
        [0, 0, 0, 0, 0, 0, 0, 0, 1, 1],
    ]
) - np.eye(10)

s = hafnian_sample_graph(matrix, 6, samples=10)
print(s)

Can you provide the entire error traceback when you try to run this? Don't leave anything out 😄

minseok1999 commented 1 year ago

@minseok1999 Let's stick with the simpler example here:

import numpy as np
from thewalrus.samples import hafnian_sample_graph

matrix = np.array(
    [
        [1, 1, 1, 1, 1, 1, 0, 0, 0, 0],
        [1, 1, 1, 1, 1, 1, 0, 0, 0, 0],
        [1, 1, 1, 1, 1, 1, 0, 0, 0, 0],
        [1, 1, 1, 1, 1, 1, 0, 0, 1, 0],
        [1, 1, 1, 1, 1, 1, 0, 1, 0, 0],
        [1, 1, 1, 1, 1, 1, 1, 0, 0, 0],
        [0, 0, 0, 0, 0, 1, 1, 0, 0, 0],
        [0, 0, 0, 0, 1, 0, 0, 1, 0, 0],
        [0, 0, 0, 1, 0, 0, 0, 0, 1, 1],
        [0, 0, 0, 0, 0, 0, 0, 0, 1, 1],
    ]
) - np.eye(10)

s = hafnian_sample_graph(matrix, 6, samples=10)
print(s)

Can you provide the entire error traceback when you try to run this? Don't leave anything out 😄

I have tried your simpler form and I got

OMP: Info #276: omp_set_nested routine deprecated, please use omp_set_max_active_levels instead. [[0 0 0 0 0 0 0 0 0 0] [1 0 1 1 2 0 0 1 0 0] [0 0 0 0 0 0 0 0 0 0] [0 0 0 0 0 1 1 0 0 0] [3 2 3 0 3 1 0 0 0 0] [4 4 2 1 5 4 0 0 0 0] [3 2 1 3 0 4 1 0 0 0] [1 1 1 1 1 1 0 0 0 0] [2 1 0 1 2 0 0 0 0 0] [1 0 1 0 0 0 0 0 0 0]]

without the division by zero error anymore. Why is this like this? Should I not be fidgeting with the fixed library? Only change I made to the library was to change the number of samples explicitly in the definition of hafnian_sample_graph itself 😅 And what is going on with this nested routine deprecation in this case?

isaacdevlugt commented 1 year ago

He @minseok1999, I asked some other folks internally to try and replicate the behaviour you're getting, and nobody is able to — everyone can run this code without error. A couple things you can try:

Let me know if either of those work!

minseok1999 commented 1 year ago

I have tried google colab as you advised and this works cleanly without error. And by creating a new file without redefining library myself and reinstalling thewalrus library in that new file (pip install thewalrus) I could see that it works fine and I could get result of hafnian_sample_graph without error. (Before the renewing of the library using pip install thewalrus this nested routine deprecation message persisted)

I wish I could figure out why this message "OMP: Info https://github.com/XanaduAI/thewalrus/pull/276: omp_set_nested routine deprecated, please use omp_set_max_active_levels instead" or " divide by zero error" occurs or locate where the clash happens when I copy and paste from the original library and slightly modify the definition itself in my own Jupyter Notebook, but I ain't an expert on how the library structure works in Jupyter Notebook so I'd better be satisfied with the clean result from sticking to your simpler code with renewed installation of thewalrus library😅

isaacdevlugt commented 1 year ago

Awesome! Glad you got this solved 😄. Sometimes virtual environments get "messed up" and need a refresh just like anything else would 😅.