diff --git a/thewalrus/decompositions.py b/thewalrus/decompositions.py index 43967077..ad0e767b 100644 --- a/thewalrus/decompositions.py +++ b/thewalrus/decompositions.py @@ -210,6 +210,7 @@ def takagi(A, svd_order=True): list_vals.sort(reverse=svd_order) sorted_ls, permutation = zip(*list_vals) return np.array(sorted_ls), Uc[:, np.array(permutation)] + # Find the element with the largest absolute value pos = np.unravel_index(np.argmax(np.abs(A)), (n, n)) # Use it to find whether the input is a global phase times a real matrix @@ -223,7 +224,7 @@ def takagi(A, svd_order=True): U = u @ sqrtm((v @ np.conjugate(u)).T) # The line above could be simplifed to the line below if the product v @ np.conjugate(u) is diagonal # Which it should be according to Caves http://info.phys.unm.edu/~caves/courses/qinfo-s17/lectures/polarsingularAutonne.pdf - # U = u * np.sqrt(0j + np.diag(np.conjugate(u) @ v)) + # U = u * np.sqrt(0j + np.diag(v @ np.conjugate(u))) # This however breaks test_degenerate if svd_order is False: return d[::-1], U[:, ::-1] diff --git a/thewalrus/symplectic.py b/thewalrus/symplectic.py index 2583103c..0321e518 100644 --- a/thewalrus/symplectic.py +++ b/thewalrus/symplectic.py @@ -438,39 +438,6 @@ def is_symplectic(S, rtol=1e-05, atol=1e-08): return np.allclose(S.T @ Omega @ S, Omega, rtol=rtol, atol=atol) -def autonne(A, rtol=1e-05, atol=1e-08, svd_order=True): - r"""Autonne-Takagi decomposition of a complex symmetric (not Hermitian!) matrix. - - Args: - A (array): square, symmetric matrix - rtol (float): the relative tolerance parameter between ``A`` and ``A.T`` - atol (float): the absolute tolerance parameter between ``A`` and ``A.T`` - svd_order (boolean): whether to return result by ordering the singular values of ``A`` in descending (``True``) or asceding (``False``) order. - - Returns: - tuple[array, array]: (r, U), where r are the singular values, - and U is the Autonne-Takagi unitary, such that :math:`A = U \diag(r) U^T`. - """ - n, m = A.shape - if n != m: - raise ValueError("The input matrix is not square") - if not np.allclose(A, A.T, rtol=rtol, atol=atol): - raise ValueError("The input matrix is not symmetric") - Areal = A.real - Aimag = A.imag - - B = np.empty((2 * n, 2 * n)) - B[:n, :n] = Areal - B[n : 2 * n, :n] = Aimag - B[:n, n : 2 * n] = Aimag - B[n : 2 * n, n : 2 * n] = -Areal - vals, vects = np.linalg.eigh(B) - U = vects[:n, n : 2 * n] + 1j * vects[n : 2 * n, n : 2 * n] - if svd_order: - return (vals[n : 2 * n])[::-1], U[:, ::-1] - return vals[n : 2 * n], U - - def xxpp_to_xpxp(S): """Permutes the entries of the input from xxpp ordering to xpxp ordering. diff --git a/thewalrus/tests/test_decompositions.py b/thewalrus/tests/test_decompositions.py index 53975dd2..ab7deeff 100644 --- a/thewalrus/tests/test_decompositions.py +++ b/thewalrus/tests/test_decompositions.py @@ -396,6 +396,25 @@ def test_real_degenerate(): assert np.allclose(U @ np.diag(rl) @ U.T, mat) +@pytest.mark.parametrize("n", [5, 10, 50]) +@pytest.mark.parametrize("datatype", [np.complex128, np.float64]) +@pytest.mark.parametrize("svd_order", [True, False]) +def test_autonne_takagi(n, datatype, svd_order): + """Checks the correctness of the Autonne decomposition function""" + if datatype is np.complex128: + A = np.random.rand(n, n) + 1j * np.random.rand(n, n) + if datatype is np.float64: + A = np.random.rand(n, n) + A += A.T + r, U = takagi(A, svd_order=svd_order) + assert np.allclose(A, U @ np.diag(r) @ U.T) + assert np.all(r >= 0) + if svd_order is True: + assert np.all(np.diff(r) <= 0) + else: + assert np.all(np.diff(r) >= 0) + + @pytest.mark.parametrize("size", [10, 20, 100]) def test_flat_phase(size): """Test that the correct decomposition is obtained even if the first entry is 0""" diff --git a/thewalrus/tests/test_symplectic.py b/thewalrus/tests/test_symplectic.py index 6afb0edd..40ff800b 100644 --- a/thewalrus/tests/test_symplectic.py +++ b/thewalrus/tests/test_symplectic.py @@ -12,7 +12,7 @@ # See the License for the specific language governing permissions and # limitations under the License. """Symplectic tests""" -# pylint: disable=no-self-use, assignment-from-no-return +# pylint: disable=assignment-from-no-return import pytest import numpy as np @@ -850,39 +850,6 @@ def test_sympmat(n): assert np.all(X == res) -@pytest.mark.parametrize("n", [5, 10, 50]) -@pytest.mark.parametrize("datatype", [np.complex128, np.float64]) -@pytest.mark.parametrize("svd_order", [True, False]) -def test_autonne(n, datatype, svd_order): - """Checks the correctness of the Autonne decomposition function""" - if datatype is np.complex128: - A = np.random.rand(n, n) + 1j * np.random.rand(n, n) - if datatype is np.float64: - A = np.random.rand(n, n) - A += A.T - r, U = symplectic.autonne(A, svd_order=svd_order) - assert np.allclose(A, U @ np.diag(r) @ U.T) - assert np.all(r >= 0) - if svd_order is True: - assert np.all(np.diff(r) <= 0) - else: - assert np.all(np.diff(r) >= 0) - - -def test_autonne_error(): - """Tests the value errors of autonne""" - n = 10 - m = 20 - A = np.random.rand(n, m) - with pytest.raises(ValueError, match="The input matrix is not square"): - symplectic.autonne(A) - n = 10 - m = 10 - A = np.random.rand(n, m) - with pytest.raises(ValueError, match="The input matrix is not symmetric"): - symplectic.autonne(A) - - class TestPhaseSpaceFunctions: """Tests for the shared phase space operations"""