diff --git a/thewalrus/decompositions.py b/thewalrus/decompositions.py index c08d6d67..72b7aef0 100644 --- a/thewalrus/decompositions.py +++ b/thewalrus/decompositions.py @@ -125,9 +125,9 @@ def blochmessiah(S): S (array[float]): 2N x 2N real symplectic matrix Returns: - tupple(array[float], : orthogonal symplectic matrix uff - array[float], : diagional matrix dff - array[float]) : orthogonal symplectic matrix vff + tuple(array[float], : orthogonal symplectic matrix uff + array[float], : diagonal matrix dff + array[float]) : orthogonal symplectic matrix vff """ N, _ = S.shape @@ -148,16 +148,11 @@ def blochmessiah(S): alpha = Unitary[0 : N // 2, 0 : N // 2] beta = Sig[0 : N // 2, N // 2 : N] # Bloch-Messiah in this Basis - u2, d2, v2 = np.linalg.svd(beta) + d2, takagibeta = takagi(beta) sval = np.arcsinh(d2) - takagibeta = u2 @ sqrtm(np.conjugate(u2).T @ (v2.T)) - uf = np.block([[takagibeta, 0 * takagibeta], [0 * takagibeta, np.conjugate(takagibeta)]]) - vf = np.block( - [ - [np.conjugate(takagibeta).T @ alpha, 0 * takagibeta], - [0 * takagibeta, np.conjugate(np.conjugate(takagibeta).T @ alpha)], - ] - ) + uf = block_diag(takagibeta, takagibeta.conj()) + blc = np.conjugate(takagibeta).T @ alpha + vf = block_diag(blc, blc.conj()) df = np.block( [ [np.diag(np.cosh(sval)), np.diag(np.sinh(sval))], @@ -172,3 +167,67 @@ def blochmessiah(S): vff = np.real_if_close(vff) uff = np.real_if_close(uff) return uff, dff, vff + + +def takagi(A, svd_order=True): + r"""Autonne-Takagi decomposition of a complex symmetric (not Hermitian!) matrix. + Note that the input matrix is internally symmetrized. If the input matrix is indeed symmetric this leaves it unchanged. + See `Carl Caves note. `_ + + Args: + A (array): square, symmetric matrix + svd_order (boolean): whether to return result by ordering the singular values of ``A`` in descending (``True``) or ascending (``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") + # Here we force symmetrize the matrix + A = 0.5 * (A + A.T) + + A = np.real_if_close(A) + + if np.allclose(A, 0): + return np.zeros(n), np.eye(n) + + if np.isrealobj(A): + # If the matrix A is real one can be more clever and use its eigendecomposition + ls, U = np.linalg.eigh(A) + U = U / np.exp(1j * np.angle(U)[0]) + vals = np.abs(ls) # These are the Takagi eigenvalues + phases = -np.ones(vals.shape[0], dtype=np.complex128) + for j, l in enumerate(ls): + if np.allclose(l, 0) or l > 0: + phases[j] = 1 + phases = np.sqrt(phases) + Uc = U @ np.diag(phases) # One needs to readjust the phases + signs = np.sign(Uc.real)[0] + for k, s in enumerate(signs): + if np.allclose(s, 0): + signs[k] = 1 + Uc = np.real_if_close(Uc / signs) + list_vals = [(vals[i], i) for i in range(len(vals))] + # And also rearrange the unitary and values so that they are decreasingly ordered + list_vals.sort(reverse=svd_order) + sorted_ls, permutation = zip(*list_vals) + return np.array(sorted_ls), Uc[:, np.array(permutation)] + + phi = np.angle(A[0, 0]) + Amr = np.real_if_close(np.exp(-1j * phi) * A) + if np.isrealobj(Amr): + vals, U = takagi(Amr, svd_order=svd_order) + return vals, U * np.exp(1j * phi / 2) + + u, d, v = np.linalg.svd(A) + 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(v @ np.conjugate(u))) + # This however breaks test_degenerate + if svd_order is False: + return d[::-1], U[:, ::-1] + return d, U diff --git a/thewalrus/tests/test_decompositions.py b/thewalrus/tests/test_decompositions.py index 0161be7a..630eb32a 100644 --- a/thewalrus/tests/test_decompositions.py +++ b/thewalrus/tests/test_decompositions.py @@ -1,4 +1,4 @@ -# Copyright 2019 Xanadu Quantum Technologies Inc. +# Copyright 2023 Xanadu Quantum Technologies Inc. # Licensed under the Apache License, Version 2.0 (the "License"); # you may not use this file except in compliance with the License. @@ -19,7 +19,7 @@ from thewalrus.random import random_interferometer as haar_measure from thewalrus.random import random_symplectic -from thewalrus.decompositions import williamson, blochmessiah +from thewalrus.decompositions import williamson, blochmessiah, takagi from thewalrus.symplectic import sympmat as omega from thewalrus.quantum.gaussian_checks import is_symplectic @@ -253,3 +253,144 @@ def test_transform(self, passive, create_transform, tol): assert np.allclose(O1.T @ O @ O1, O, atol=tol, rtol=0) assert np.allclose(O2.T @ O @ O2, O, atol=tol, rtol=0) assert np.allclose(S @ O @ S.T, O, atol=tol, rtol=0) + + +@pytest.mark.parametrize("n", [5, 10, 50]) +@pytest.mark.parametrize("datatype", [np.complex128, np.float64]) +@pytest.mark.parametrize("svd_order", [True, False]) +def test_takagi(n, datatype, svd_order): + """Checks the correctness of the Takagi 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("n", [5, 10, 50]) +@pytest.mark.parametrize("datatype", [np.complex128, np.float64]) +@pytest.mark.parametrize("svd_order", [True, False]) +@pytest.mark.parametrize("half_rank", [0, 1]) +@pytest.mark.parametrize("phase", [0, 1]) +def test_degenerate(n, datatype, svd_order, half_rank, phase): + """Tests Takagi produces the correct result for very degenerate cases""" + nhalf = n // 2 + diags = [half_rank * np.random.rand()] * nhalf + [np.random.rand()] * (n - nhalf) + if datatype is np.complex128: + U = haar_measure(n) + if datatype is np.float64: + U = np.exp(1j * phase) * haar_measure(n, real=True) + A = U @ np.diag(diags) @ U.T + r, U = takagi(A, svd_order=svd_order) + assert np.allclose(A, U @ np.diag(r) @ U.T) + assert np.allclose(U @ U.T.conj(), np.eye(n)) + 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_zeros(): + """Verify that the Takagi decomposition returns a zero vector and identity matrix when + input a matrix of zeros""" + dim = 4 + a = np.zeros((dim, dim)) + rl, U = takagi(a) + assert np.allclose(rl, np.zeros(dim)) + assert np.allclose(U, np.eye(dim)) + + +def test_takagi_error(): + """Tests the value errors of Takagi""" + n = 10 + m = 11 + A = np.random.rand(n, m) + with pytest.raises(ValueError, match="The input matrix is not square"): + takagi(A) + + +def test_real_degenerate(): + """Verify that the Takagi decomposition returns a matrix that is unitary and results in a + correct decomposition when input a real but highly degenerate matrix. This test uses the + adjacency matrix of a balanced tree graph.""" + + vals = [ + 1, + 2, + 31, + 34, + 35, + 62, + 67, + 68, + 94, + 100, + 101, + 125, + 133, + 134, + 157, + 166, + 167, + 188, + 199, + 200, + 220, + 232, + 233, + 251, + 265, + 266, + 283, + 298, + 299, + 314, + 331, + 332, + 346, + 364, + 365, + 377, + 397, + 398, + 409, + 430, + 431, + 440, + 463, + 464, + 472, + 503, + 535, + 566, + 598, + 629, + 661, + 692, + 724, + 755, + 787, + 818, + 850, + 881, + 913, + 944, + ] + mat = np.zeros([31 * 31]) + mat[vals] = 1 + mat = mat.reshape(31, 31) + # The lines above are equivalent to: + # import networkx as nx + # g = nx.balanced_tree(2, 4) + # a = nx.to_numpy_array(g) + rl, U = takagi(mat) + assert np.allclose(U @ U.conj().T, np.eye(len(mat))) + assert np.allclose(U @ np.diag(rl) @ U.T, mat)