Default method fails on poorly conditioned matrices

[jacobhilton]

Before posting a bug report

• [x] I have searched exisisting GitHub issues to make sure the issue does not already exist.

Expected behavior

Calling hafnian on a poorly conditioned matrix should produce something close to the correct answer.

Actual behavior

An incorrect answer is produced. For example, on an 8 x 8 matrix of ones with a single off-diagonal entry equal to 1,000,000, an answer of 0 is produced instead of the correct answer of 15,000,090.

Reproduces how often

The behavior seems to be deterministic.

System information

Python version:            3.10.8
Platform info:             macOS-13.4.1-arm64-arm-64bit
Installation path:         /Users/jacob/miniconda3/lib/python3.10/site-packages/thewalrus
The Walrus version:        0.20.0
Numpy version:             1.24.2
Scipy version:             1.10.1
SymPy version:             1.11.1
Numba version:             0.57.0

Source code

import numpy as np
from thewalrus import hafnian

mat = np.ones((8, 8))
mat[0, 1] = mat[1, 0] = 1e6
print(hafnian(mat))  # Produces wrong answer
print(hafnian(mat, method="inclexcl"))  # Produces different but also wrong answer
print(hafnian(mat, method="recursive"))  # Works fine

Tracebacks

No response

Additional information

I suspect that the trace algorithm is unstable for poorly conditioned matrices. Perhaps it would be good to add a check to revert to the recursive algorithm in this case.

[isaacdevlugt]

Hey @jacobhilton! We're taking a look at this — will get back to you as soon as we can!

[isaacdevlugt]

@jacobhilton we think you may have found a bug, but it appears to be a relatively extreme edge case. There is likely a deeper numerical research project that needs to be done here.

It might be worth you checking out this page on matrix conditioning or elsewhere to see if you can find a criterion for which this can be dealt with as a special case.

[jacobhilton]

Thanks. Upon reflection, I think the condition number is probably not a good criterion. For example, all three methods work fine on a matrix of ones, which is singular. I don't know exactly what criterion would make sense as I don't understand the trace algorithm well enough and haven't tried to debug where it's going wrong in this particular case. I wouldn't be surprised if there's some intermediate matrix involved that needs to be well conditioned. For my use case, I will probably just stick to the recursive method – I just thought you'd want to know about the issue.

[isaacdevlugt]

@jacobhilton thanks so much for bringing this to our attention! I think we'll leave the issue open, as we didn't / haven't come to a satisfactory solution. That said, we're not too concerned with the behaviour since it's an edge case.