From e8cf951d3e7e3754d80001f2873481f299dee42e Mon Sep 17 00:00:00 2001 From: Ryosuke Noro <64354442+RyosukeNORO@users.noreply.github.com> Date: Wed, 3 Jul 2024 17:05:27 -0400 Subject: [PATCH 1/5] This is a fix of a sqrtm problem --- .github/CHANGELOG.md | 2 ++ thewalrus/decompositions.py | 5 +++-- 2 files changed, 5 insertions(+), 2 deletions(-) diff --git a/.github/CHANGELOG.md b/.github/CHANGELOG.md index 66dc415a..d6b883e4 100644 --- a/.github/CHANGELOG.md +++ b/.github/CHANGELOG.md @@ -15,6 +15,8 @@ ### Bug fixes +* Added tolerance parameter to `deecompositions.takagi` to prevent from incorrect Autonne-Takagi decomposition, which happened with some matrices. + ### Documentation ### Contributors diff --git a/thewalrus/decompositions.py b/thewalrus/decompositions.py index 78f2575c..51a299fe 100644 --- a/thewalrus/decompositions.py +++ b/thewalrus/decompositions.py @@ -152,7 +152,7 @@ def blochmessiah(S): return O, D, Q -def takagi(A, svd_order=True): +def takagi(A, svd_order=True, tol=1e-16): 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. @@ -162,6 +162,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. + tol (float): tolerance parameter to kill the small values in the maatix A. This was installed because takagi decomposition did not work at some matrix A. Returns: tuple[array, array]: (r, U), where r are the singular values, @@ -174,7 +175,7 @@ def takagi(A, svd_order=True): # Here we build a Symmetric matrix from the top right triangular part A = np.triu(A) + np.triu(A, k=1).T - A = np.real_if_close(A) + A = np.where(np.abs(A.real) <= tol, 0, A.real) + np.where(np.abs(A.imag) <= tol, 0, A.imag) if np.allclose(A, 0): return np.zeros(n), np.eye(n) From 283ae0c28c6da96a3d1e6f83a0e39be010d6b4ec Mon Sep 17 00:00:00 2001 From: Ryosuke Noro <64354442+RyosukeNORO@users.noreply.github.com> Date: Thu, 4 Jul 2024 10:00:04 -0400 Subject: [PATCH 2/5] Update decompositions.py --- thewalrus/decompositions.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/thewalrus/decompositions.py b/thewalrus/decompositions.py index 51a299fe..57239250 100644 --- a/thewalrus/decompositions.py +++ b/thewalrus/decompositions.py @@ -162,7 +162,7 @@ def takagi(A, svd_order=True, tol=1e-16): 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. - tol (float): tolerance parameter to kill the small values in the maatix A. This was installed because takagi decomposition did not work at some matrix A. + tol (float): tolerance parameter to kill small values in the matrix A. This was installed because takagi decomposition did not work at some matrix A. Returns: tuple[array, array]: (r, U), where r are the singular values, @@ -175,6 +175,7 @@ def takagi(A, svd_order=True, tol=1e-16): # Here we build a Symmetric matrix from the top right triangular part A = np.triu(A) + np.triu(A, k=1).T + # Kill small values in the matrix A. This was installed because takagi decomposition did not work at some matrix A with small values. A = np.where(np.abs(A.real) <= tol, 0, A.real) + np.where(np.abs(A.imag) <= tol, 0, A.imag) if np.allclose(A, 0): From b65026ac57f3d46aa83c49f3f1244e3e423a2b46 Mon Sep 17 00:00:00 2001 From: Ryosuke Noro <64354442+RyosukeNORO@users.noreply.github.com> Date: Thu, 4 Jul 2024 10:12:44 -0400 Subject: [PATCH 3/5] Update prepare_cov.py --- thewalrus/internal_modes/prepare_cov.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/thewalrus/internal_modes/prepare_cov.py b/thewalrus/internal_modes/prepare_cov.py index fc7e12a2..b4b368ec 100644 --- a/thewalrus/internal_modes/prepare_cov.py +++ b/thewalrus/internal_modes/prepare_cov.py @@ -230,7 +230,7 @@ def state_prep(eps, W, thresh=1e-4, hbar=2): def LO_overlaps(chis, LO_shape): r""" - Computes the overlap integral between the orthonormal moes and the local oscillator shape + Computes the overlap integral between the orthonormal modes and the local oscillator shape Args: chis (list[array]): list of the temporal functions for the new orthonormal basis From 911719807a36852199ab554091d753e9e7780f18 Mon Sep 17 00:00:00 2001 From: Ryosuke Noro <64354442+RyosukeNORO@users.noreply.github.com> Date: Thu, 4 Jul 2024 10:29:52 -0400 Subject: [PATCH 4/5] added the related pull request number into CHANGELOG.md. --- .github/CHANGELOG.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/CHANGELOG.md b/.github/CHANGELOG.md index d6b883e4..09b945ba 100644 --- a/.github/CHANGELOG.md +++ b/.github/CHANGELOG.md @@ -15,7 +15,7 @@ ### Bug fixes -* Added tolerance parameter to `deecompositions.takagi` to prevent from incorrect Autonne-Takagi decomposition, which happened with some matrices. +* Added tolerance parameter to `deecompositions.takagi` to prevent from incorrect Autonne-Takagi decomposition, which happened with some matrices. [(#393)](https://github.com/XanaduAI/thewalrus/pull/393) ### Documentation From 1609f1376e8afcfe192ca7c9209aa63dcfcb4b5b Mon Sep 17 00:00:00 2001 From: Ryosuke Noro <64354442+RyosukeNORO@users.noreply.github.com> Date: Thu, 4 Jul 2024 14:24:42 -0400 Subject: [PATCH 5/5] Update thewalrus/decompositions.py Co-authored-by: Eli Bourassa <53090166+elib20@users.noreply.github.com> --- thewalrus/decompositions.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/thewalrus/decompositions.py b/thewalrus/decompositions.py index 57239250..3ca4dd29 100644 --- a/thewalrus/decompositions.py +++ b/thewalrus/decompositions.py @@ -176,7 +176,8 @@ def takagi(A, svd_order=True, tol=1e-16): A = np.triu(A) + np.triu(A, k=1).T # Kill small values in the matrix A. This was installed because takagi decomposition did not work at some matrix A with small values. - A = np.where(np.abs(A.real) <= tol, 0, A.real) + np.where(np.abs(A.imag) <= tol, 0, A.imag) + rtol = tol * np.amax(np.abs(A)) + A = np.where(np.abs(A.real) <= rtol, 0, A.real) + np.where(np.abs(A.imag) <= rtol, 0, A.imag) if np.allclose(A, 0): return np.zeros(n), np.eye(n)