Skip to content

Commit

Permalink
Merge branch 'master' into internal_modes
Browse files Browse the repository at this point in the history
  • Loading branch information
RyosukeNORO authored Jul 8, 2024
2 parents f7358e6 + 0e51518 commit 2b249e9
Show file tree
Hide file tree
Showing 3 changed files with 47 additions and 1 deletion.
2 changes: 2 additions & 0 deletions .github/CHANGELOG.md
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,8 @@

### Bug fixes

* Add the calculation method of `takagi` when the matrix is diagonal. [(#394)](https://github.com/XanaduAI/thewalrus/pull/394)

### Documentation

### Contributors
Expand Down
17 changes: 16 additions & 1 deletion thewalrus/decompositions.py
Original file line number Diff line number Diff line change
Expand Up @@ -152,7 +152,8 @@ def blochmessiah(S):
return O, D, Q


def takagi(A, svd_order=True):
def takagi(A, svd_order=True, rtol=1e-16):
# pylint: disable=too-many-return-statements
r"""Autonne-Takagi decomposition of a complex symmetric (not Hermitian!) matrix.
Note that the input matrix is internally symmetrized by taking its upper triangular part.
If the input matrix is indeed symmetric this leaves it unchanged.
Expand All @@ -162,6 +163,7 @@ def takagi(A, svd_order=True):
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.
rtol (float): the relative tolerance parameter used in ``np.allclose`` when judging if the matrix is diagonal or not. Default to 1e-16.
Returns:
tuple[array, array]: (r, U), where r are the singular values,
Expand Down Expand Up @@ -202,6 +204,19 @@ def takagi(A, svd_order=True):
vals, U = takagi(Amr, svd_order=svd_order)
return vals, U * np.exp(1j * phi / 2)

# If the matrix is diagonal, Takagi decomposition is easy
if np.allclose(A, np.diag(np.diag(A)), rtol=rtol):
d = np.diag(A)
l = np.abs(d)
idx = np.argsort(l)
d = d[idx]
l = l[idx]
U = np.diag(np.exp(1j * 0.5 * np.angle(d)))
U = U[::-1, :]
if svd_order:
return l[::-1], U[:, ::-1]
return l, U

u, d, v = np.linalg.svd(A)
U = u @ sqrtm((v @ np.conjugate(u)).T)
if svd_order is False:
Expand Down
29 changes: 29 additions & 0 deletions thewalrus/tests/test_decompositions.py
Original file line number Diff line number Diff line change
Expand Up @@ -324,6 +324,35 @@ def test_takagi_error():
takagi(A)


@pytest.mark.parametrize("svd_order", [True, False])
def test_takagi_diagonal_matrix(svd_order):
"""Test the takagi decomposition works well for a specific matrix that was not decomposed accurately in a previous implementation.
See more info in PR #393 (https://github.com/XanaduAI/thewalrus/pull/393)"""
A = np.array(
[
[
-8.4509484628125742e-01 + 1.0349426984742664e-16j,
6.3637197288239186e-17 - 7.4398922703555097e-33j,
2.6734481396039929e-32 + 1.7155650257063576e-35j,
],
[
6.3637197288239186e-17 - 7.4398922703555097e-33j,
-2.0594021562561332e-01 + 2.2863956908382538e-17j,
-5.8325863096557049e-17 + 1.6949718400585382e-18j,
],
[
2.6734481396039929e-32 + 1.7155650257063576e-35j,
-5.8325863096557049e-17 + 1.6949718400585382e-18j,
4.4171453199503476e-02 + 1.0022350742842835e-02j,
],
]
)
d, U = takagi(A, svd_order=svd_order)
assert np.allclose(A, U @ np.diag(d) @ U.T)
assert np.allclose(U @ np.conjugate(U).T, np.eye(len(U)))
assert np.all(d >= 0)


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
Expand Down

0 comments on commit 2b249e9

Please sign in to comment.