Incompatibility with Numpy 1.25

[jakelishman]

Environment

• Qiskit Terra version: 0.24.1
• Python version: 3.10
• Operating system: macOS 13.3.1

What is happening?

The test suite and CI are failing because of new deprecation warnings issued by Numpy 1.25 (removing some unnecessary aliases and deprecating some questionable array-to-scalar conversions that we first encountered in #10120). There also seems to be some sort of behavioural change that's affecting the isometry tests and making them flaky, which is harder to pin down.

How can we reproduce the issue?

Run CI with Numpy 1.25.

What should happen?

No warnings, all tests passing.

Any suggestions?

No response

[jakelishman]

I'm not 100% sure what's going on with the isometry tests - the failures seem to have questionably large absolute differences between the expected and resultant matrices. It's plausible that it's related to changes in sin and cos in Numpy 1.25: numpy/numpy#23399, but I'm still looking into it.

[jakelishman]

Here's a reproducer for flaky test behaviour with Numpy 1.25:

import numpy as np
from qiskit.quantum_info import random_unitary, Operator
from qiskit.extensions.quantum_initializer import Isometry

iso = random_unitary(16, seed=278663).data
np.testing.assert_allclose(Operator(Isometry(iso, 0, 0)), Operator(iso))

There are no ancilla qubits, so I think the unitary from the isometry is fully defined and there should be no deviations beyond standard floating-point stuff. Instead, the Operator(Isometry(iso, 0, 0)) part is showing some very nontrivial non-deterministic variance in the off-diagonal elements each time it's called.

[jakelishman]

I've narrowed at least one component of this a little bit more. This might not be the only issue, but I do know:

import numpy as np
from numpy import array
from qiskit.extensions.quantum_initializer import UCGate

np.set_printoptions(floatmode="unique")

a1 = [array([[ 0.668014861563423-0.3346551458465917j,  0.277967131510269-0.6037353326489598j],
       [-0.277967131510269-0.6037353326489598j,  0.668014861563423+0.3346551458465917j]]), array([[ 0.2166062838266952+0.8682003952928524j ,  0.4271167766350859+0.12992709700820013j],
       [-0.4271167766350859+0.12992709700820013j,  0.2166062838266952-0.8682003952928524j ]]), array([[-0.2741070097201944 +0.8117613476501147j, -0.01821474636319918+0.5153417164352673j],
       [ 0.01821474636319918+0.5153417164352673j, -0.2741070097201944 -0.8117613476501147j]]), array([[ 0.1623765358315201 +0.9716678872200901j,  0.13410238340678943-0.107293654617689j ],
       [-0.13410238340678943-0.107293654617689j ,  0.1623765358315201 -0.9716678872200901j]]), array([[ 0.8245314897717796  +0.5297803036914427j , -0.014618052150304811-0.19815893809430707j],
       [ 0.014618052150304811-0.19815893809430707j,  0.8245314897717796  -0.5297803036914427j ]]), array([[ 0.7996715550698764 +0.5389276848385285j, -0.26190637639918435-0.0385668838274234j],
       [ 0.26190637639918435-0.0385668838274234j,  0.7996715550698764 -0.5389276848385285j]]), array([[-0.6951654137555511-0.4727347508596223j , -0.4217756893779492+0.33966479166601926j],
       [ 0.4217756893779492+0.33966479166601926j, -0.6951654137555511+0.4727347508596223j ]]), array([[-0.077850612760931  +0.5820993526489561j,  0.06379232106309757+0.8068644034233984j],
       [-0.06379232106309757+0.8068644034233984j, -0.077850612760931  -0.5820993526489561j]])]
uc = UCGate(a1, up_to_diagonal=True)

# As of Numpy 1.25, this call is no longer bit-for-bit deterministic:
uc._get_diagonal()

That call at the end used to produce deterministic output. Now, the various values can differ between calls by an atol of somewhere in the region of 2e-16. This then threads itself through the rest of the functions, occasionally flipping some global phases, and producing different matrices. There's a call to np.linalg.eig in the middle of that above call, which I'm always suspicious of, but the 1.25 changelog doesn't mention anything obvious about how that might be affected.

This doesn't explain to me, however, why the operator for the output isometry differs significantly from the input. Within a single construction, I don't know why minor differences in the FLOPs of the calculation can produce such different output.

[jakelishman]

I don't know if this is the root cause, but I've tracked down a seeming flakiness in the multiply ufunc with complex numbers in Numpy 1.25. On my machine, the file

import sys
import numpy

a = [
    numpy.array([
        [0.37407431362525206 + 0.7981674467868007j, 0.4497329415302707 + 0.1440049168541349j],
        [0.44973294153027066 - 0.14400491685413508j, -0.37407431362525156 + 0.7981674467868008j],
    ])
    for _ in [None] * 2
]
b = 0.7071067811865476 - 0.7071067811865475j

if (a[0] * b)[1, 1] - (a[1] * b)[1, 1]:
    sys.exit(1)

intermittently exits in the failure condition, with the imaginary components differing by 1ULP. I couldn't reproduce that behaviour with Numpy 1.24.3 (or 1.25.0 without AVX2).

I don't really understand how that behaviour could cause a correctness issue in the Isometry decomposition itself - as far as I can tell it shouldn't attempt to recalculate the exact same FLOPs at any other time, so the decomposition should come out maybe differently each time, but correct - but it at least seems to be related.

I've filed the above issue as numpy/numpy#24000. I'm not 100% sure it's a bug in Numpy, but I've chased this about as far as I can think right now.

[jakelishman]

This comment is meant to serve as a documentation point for workarounds for UCGate/multiplixer-gate synthesis problems with NumPy 1.25 and 1.26 on Intel Macs with AVX2 CPU extensions enabled. The comments above this one describe a problem with these versions of NumPy in these conditions, which causes severe problems during our synthesis of multiplexed gates, most commonly encountered in Qiskit by calling UnitaryGate.control().

There are a few options for working around the bug, in order from easiest to most complex:

• Install an unaffected version of NumPy. At this point in time, we know that of PyPI releases: the 1.24 series and everything before it are good; the entire 1.25 series is bad; the current most-recent release 1.26.0 is bad. I haven't checked the versions available on conda or conda-forge. You can track upstream progress on fixing the bug in BUG: Multiply complex128 ufunc inconsistent with AVX2 numpy/numpy#24000.

• Disable NumPy's use of the AVX2 CPU extensions at runtime. To do this, set the environment variable NPY_DISABLE_CPU_FEATURES=AVX2 before importing NumPy. Typically this would mean launching python with that environment variable set, such as by running NPY_DISABLE_CPU_FEATURES=AVX2 python.

• Install a version of the affected NumPy releases that was compiled using a compiler that does not exhibit the AVX2 bug. See BUG: Multiply complex128 ufunc inconsistent with AVX2 numpy/numpy#24000 for details on the necessary compilers.

  If Issue warning on known-bad combinations for UCGate synthesis #10906 has merged, and if you are sure that your NumPy installation is good, you can set the environment variable QISKIT_CMUL_AVX2_GOOD_NUMPY=1 before importing Qiskit to disable runtime warnings about the bad version.