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

qiskit/quantum_info/synthesis/two_qubit_decompose.py

@@ -169,7 +169,7 @@ def __init__(self, unitary_matrix, eps=1e-15):
         for i in range(100):  # FIXME: this randomized algorithm is horrendous
             state = np.random.default_rng(i)
             M2real = state.normal()*M2.real + state.normal()*M2.imag
-            _, P = la.eigh(M2real)
+            _, P = np.linalg.eigh(M2real)
             D = P.T.dot(M2).dot(P).diagonal()
             if np.allclose(P.dot(np.diag(D)).dot(P.T), M2, rtol=1.0e-13, atol=1.0e-13):
                 break