Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Round M2 dot product in Weyl decomposition #4835

Merged
merged 8 commits into from
Aug 3, 2020

Conversation

mtreinish
Copy link
Member

Summary

Differences in floating point precision between different systems and
environments was causing inconsistent decompositions for decompositions
with 0s in the unitary would have differing results because of very
subtle floating point precision errors in the results. For example one
system would return 0 for a value and another would return
-1.55582133e-19 (which is essentially 0). This would have compound
effects in later stages because the signs were different. This commit
attempts to fix cases like this by rounding to 13 decimal places to
prevent issues like this in the future.

Details and comments

@mtreinish mtreinish requested a review from a team as a code owner July 30, 2020 21:26
Differences in floating point precision between different systems and
environments was causing inconsistent decompositions for decompositions
with 0s in the unitary would have differing results because of very
subtle floating point precision errors in the results. For example one
system would return 0 for a value and another would return
-1.55582133e-19 (which is essentially 0). This would have compound
effects in later stages because the signs were different. This commit
attempts to fix cases like this by rounding to 13 decimal places to
prevent issues like this in the future.
@mtreinish mtreinish added this to the 0.15 milestone Jul 31, 2020
@levbishop
Copy link
Member

In general I'm loth to add these kind of roundings because it usually indicates an underlying problem being papered over. At the very least when I come across code with round-to-epsilons in it, that indicates to me scars from previous battles with numeric bugs and makes me very reluctant to make changes in order to avoid reopening old wounds. To reduce such fears can we at least capture the problematic cases as a test instance and add a comment with a pointer at the lines with the rounding?

I've been playing around off and on with a couple other implementations for Weyl decomposition to eliminate the ugly randomization anyway, so with a test case I'll be able to check if any of them is able to solve address this issue without rounding.

@mtreinish
Copy link
Member Author

In general I'm loth to add these kind of roundings because it usually indicates an underlying problem being papered over. At the very least when I come across code with round-to-epsilons in it, that indicates to me scars from previous battles with numeric bugs and makes me very reluctant to make changes in order to avoid reopening old wounds. To reduce such fears can we at least capture the problematic cases as a test instance and add a comment with a pointer at the lines with the rounding?

I've been playing around off and on with a couple other implementations for Weyl decomposition to eliminate the ugly randomization anyway, so with a test case I'll be able to check if any of them is able to solve address this issue without rounding.

Sure I'll add a test case for the example that prompted me to add this. I definitely know that pain of debugging numeric bugs. I spent way too much time digging through things here to figure out that it was the CX gate decomposition in the basis being subtly off of zero here for some envs that was causing different decompositions on different systems and not an issue with the random unitary being tested..

@nonhermitian
Copy link
Contributor

Not sure what the concern is here. Setting was is a floating point zero to exactly zero is fairly common. Especially if the sign is important later.

@levbishop
Copy link
Member

levbishop commented Jul 31, 2020

Agree that it's common. IMO it should be much less common and used as a last resort.
It increases the error of the routine, by in this case factor of around 1E6. If the results of that routine feed into other routines with similar epsilon tolerances then those will have to broaden their tolerances. In the worst case you feed the output of this routine into a downstream routine that also needs to round-to-epsilon and because of this needs to broaden its tolerances correspondingly.
A factor of 1E6 here, a factor of 1E6 there, pretty soon you're talking real money ;-)

This commit adds a test for a cx decomposition and checks that
approxomite 0s are not resulting in a negative number which could result
in inconsistent results. In the process of adding that test a new place
was found where precision issues could result in an inconsistency, this
is fixed too.
@nonhermitian nonhermitian self-requested a review August 2, 2020 10:57
Copy link
Member

@ajavadia ajavadia left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think it's important to make the output reproducible on different platforms, so I'm good with this. A more robust Weyl decomposition would definitely be good.

@mergify mergify bot merged commit 8e4a366 into Qiskit:master Aug 3, 2020
@mtreinish mtreinish deleted the reduce-fp-precision branch August 3, 2020 09:43
mtreinish added a commit to mtreinish/qiskit-core that referenced this pull request Aug 4, 2020
In Qiskit#4656 we added rounding to the output of the decomp0 method to handle
a case where differing FP precision on windows environments was causing
an expected result in running two_qubit_cnot_decompose on np.eye(4)
with numpy 1.19.x installed leading to a hard failure in the qasm tests.
This seemed to reliably unblock testing and make unit tests work
reliably. However, that original fix from Qiskit#4656 was superseded by Qiskit#4835
which was a fix for a more general issue with the reproducibility of the
decompositions and reverted. Since Qiskit#4835 has merged we've been seeing an
uptick in the failure rate on the same unitary qasm test that Qiskit#4656
fixed, so the change in Qiskit#4835 was not actually sufficient for the
windows case. This commit restores the fix from Qiskit#4656 to unblock CI and
fix the reproducability of the decompositions across systems.

Fixes Qiskit#4856
kdk pushed a commit that referenced this pull request Aug 4, 2020
In #4656 we added rounding to the output of the decomp0 method to handle
a case where differing FP precision on windows environments was causing
an expected result in running two_qubit_cnot_decompose on np.eye(4)
with numpy 1.19.x installed leading to a hard failure in the qasm tests.
This seemed to reliably unblock testing and make unit tests work
reliably. However, that original fix from #4656 was superseded by #4835
which was a fix for a more general issue with the reproducibility of the
decompositions and reverted. Since #4835 has merged we've been seeing an
uptick in the failure rate on the same unitary qasm test that #4656
fixed, so the change in #4835 was not actually sufficient for the
windows case. This commit restores the fix from #4656 to unblock CI and
fix the reproducability of the decompositions across systems.

Fixes #4856
mtreinish added a commit to mtreinish/qiskit-core that referenced this pull request Oct 16, 2020
This commit removes a test Qiskit#4835 which was added to validate that the
internal values for the decomposer didn't sign flip in the specific case
that the PR was fixing. However with the use of the numpy eigh()
function these internal values no longer the same. So instead of
updating the test this just deletes it since the original case it was
catching is no longer applicable.
ajavadia pushed a commit that referenced this pull request Oct 20, 2020
* Use eigh from numpy for TwoQubitWeylDecomposition

The hermitian optimized eigen solver function in scipy has shown
to possibly produce varied results in different environments. For
example, the matrix returned by the eigh() call in the
TwoQubitWeylDecomposition class's __init__ method when called for
a CXGate:

TwoQubitWeylDecomposition(Operator(CXGate()))

was observed to be

[[ 0.70710678  0.         -0.70710678  0.        ]
 [ 0.         -0.70710678  0.          0.70710678]
 [ 0.          0.70710678  0.          0.70710678]
 [ 0.70710678  0.          0.70710678  0.        ]]

on one system/environment but

[[ 0.70710678  0.         -0.70710678  0.        ]
 [ 0.          0.70710678  0.          0.70710678]
 [ 0.         -0.70710678  0.          0.70710678]
 [ 0.70710678  0.          0.70710678  0.        ]]

on a different environment. This difference was resulting in
inconsistent decompositions being generated and made reproducing results
from the transpiler difficult. In testing the numpy equivalent
function [2] has proven to be more reliable across different
environments. This commit switches the eigh() usage in the
TwoQubitWeylDecomposition class to use the numpy function instead
this should fix the inconsistent decompositions we were seeing between
environments.

[1] https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.eigh.html
[2] https://numpy.org/doc/stable/reference/generated/numpy.linalg.eigh.html

* Remove unecessary test

This commit removes a test #4835 which was added to validate that the
internal values for the decomposer didn't sign flip in the specific case
that the PR was fixing. However with the use of the numpy eigh()
function these internal values no longer the same. So instead of
updating the test this just deletes it since the original case it was
catching is no longer applicable.

* Switch back to scipy.eigh and use evd driver

This commit reverts back to using the scipy eigh() function but sets a
the lapack driver to be fixed as evd. The numpy routine was using the
evd driver while scipy would default to use ev as the driver. Switching
the scipy driver from ev to evd seems to be the key difference between
the routines which was causing the inconsistent results with scipy.
Approaching it this way means we avoid need to bump the minimum numpy
version in the requirements list since the numpy eigh() function was
added in a fairly recent version.

* Fix lint
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants