Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Adds extra degenerate tests with nonvac nullspace #377

Merged
merged 9 commits into from
Jan 4, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 3 additions & 0 deletions .github/CHANGELOG.md
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@
* Adds the Takagi decomposition [(#363)](https://github.com/XanaduAI/thewalrus/pull/363)

* Adds the Montrealer and Loop Montrealer functions [(#363)](https://github.com/XanaduAI/thewalrus/pull/374).

### Breaking changes

### Improvements
Expand All @@ -18,6 +19,8 @@

* Improves the handling of an edge case in Takagi [(#373)](https://github.com/XanaduAI/thewalrus/pull/373).

* Adds extra tests for the Takagi decomposition [(#377)](https://github.com/XanaduAI/thewalrus/pull/377)

### Bug fixes

### Documentation
Expand Down
37 changes: 13 additions & 24 deletions thewalrus/decompositions.py
Original file line number Diff line number Diff line change
Expand Up @@ -166,7 +166,8 @@ def blochmessiah(S):

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.
Note that the input matrix is internally symmetrized by taking its upper triangular part.
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is also what np.linalg.eig and np.linalg.eigvalsh do...

If the input matrix is indeed symmetric this leaves it unchanged.
See `Carl Caves note. <http://info.phys.unm.edu/~caves/courses/qinfo-s17/lectures/polarsingularAutonne.pdf>`_

Args:
Expand All @@ -181,8 +182,8 @@ def takagi(A, svd_order=True):
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)
# 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)

Expand All @@ -192,24 +193,16 @@ def takagi(A, svd_order=True):
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)]
signs = (-1) ** (1 + np.heaviside(ls, 1))
phases = np.sqrt(np.complex128(signs))
Uc = U * phases # One needs to readjust the phases
# Find the permutation to sort in decreasing order
perm = np.argsort(vals)
# if svd_order reverse it
if svd_order:
perm = perm[::-1]
return vals[perm], Uc[:, perm]

# Find the element with the largest absolute value
pos = np.unravel_index(np.argmax(np.abs(A)), (n, n))
Expand All @@ -222,10 +215,6 @@ def takagi(A, svd_order=True):

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
17 changes: 12 additions & 5 deletions thewalrus/tests/test_decompositions.py
Original file line number Diff line number Diff line change
Expand Up @@ -274,23 +274,30 @@ def test_takagi(n, datatype, svd_order):
assert np.all(np.diff(r) >= 0)


# pylint: disable=too-many-arguments
@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):
@pytest.mark.parametrize("null_space", [0, 5, 10])
@pytest.mark.parametrize("offset", [0, 0.5])
def test_degenerate(n, datatype, svd_order, half_rank, phase, null_space, offset):
"""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)
diags = (
[half_rank * np.random.rand()] * nhalf
+ [np.random.rand() - offset] * (n - nhalf)
+ [0] * null_space
)
if datatype is np.complex128:
U = haar_measure(n)
U = haar_measure(n + null_space)
if datatype is np.float64:
U = np.exp(1j * phase) * haar_measure(n, real=True)
U = np.exp(1j * phase) * haar_measure(n + null_space, 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.allclose(U @ U.T.conj(), np.eye(n + null_space))
assert np.all(r >= 0)
if svd_order is True:
assert np.all(np.diff(r) <= 0)
Expand Down
Loading