diff --git a/.github/CHANGELOG.md b/.github/CHANGELOG.md
index ca3dce6b7..0c1bdd447 100644
--- a/.github/CHANGELOG.md
+++ b/.github/CHANGELOG.md
@@ -10,6 +10,8 @@
instead accessible via the `sf.program_utils.validate_gate_parameters` function.
[(#720)](https://github.com/XanaduAI/strawberryfields/pull/720)
+* Removes redundant decompositions already available in `thewalrus`
+
Bug fixes
Documentation
@@ -18,7 +20,7 @@
This release contains contributions from (in alphabetical order):
-Theodor Isacsson
+Theodor Isacsson, Nicolas Quesada
# Release 0.23.0 (current release)
diff --git a/requirements.txt b/requirements.txt
index bcd80965d..9dafbdc28 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -9,7 +9,7 @@ scipy==1.8.0
sympy==1.10
tensorflow==2.8.0
tensorboard==2.8.0
-thewalrus==0.19.0
+git+https://github.com/xanaduai/thewalrus.git@master
toml==0.10.2
typing-extensions==4.1.1
urllib3==1.26.8
diff --git a/strawberryfields/decompositions.py b/strawberryfields/decompositions.py
index f4f3d340d..41da27c03 100644
--- a/strawberryfields/decompositions.py
+++ b/strawberryfields/decompositions.py
@@ -16,90 +16,12 @@
used to perform gate decompositions.
"""
-from itertools import groupby
from collections import defaultdict
import numpy as np
-from scipy.linalg import block_diag, sqrtm, polar, schur
from thewalrus.quantum import adj_scaling
-from thewalrus.symplectic import sympmat, xpxp_to_xxpp
-
-
-def takagi(N, tol=1e-13, rounding=13):
- r"""Autonne-Takagi decomposition of a complex symmetric (not Hermitian!) matrix.
-
- Note that singular values of N are considered equal if they are equal after np.round(values, tol).
-
- See :cite:`cariolaro2016` and references therein for a derivation.
-
- Args:
- N (array[complex]): square, symmetric matrix N
- rounding (int): the number of decimal places to use when rounding the singular values of N
- tol (float): the tolerance used when checking if the input matrix is symmetric: :math:`|N-N^T| <` tol
-
- Returns:
- tuple[array, array]: (rl, U), where rl are the (rounded) singular values,
- and U is the Takagi unitary, such that :math:`N = U \diag(rl) U^T`.
- """
- (n, m) = N.shape
- if n != m:
- raise ValueError("The input matrix must be square")
- if np.linalg.norm(N - np.transpose(N)) >= tol:
- raise ValueError("The input matrix is not symmetric")
-
- N = np.real_if_close(N)
-
- if np.allclose(N, 0):
- return np.zeros(n), np.eye(n)
-
- if np.isrealobj(N):
- # If the matrix N is real one can be more clever and use its eigendecomposition
- l, U = np.linalg.eigh(N)
- vals = np.abs(l) # These are the Takagi eigenvalues
- phases = np.sqrt(np.complex128([1 if i > 0 else -1 for i in l]))
- Uc = U @ np.diag(phases) # One needs to readjust the phases
- list_vals = [(vals[i], i) for i in range(len(vals))]
- list_vals.sort(reverse=True)
- sorted_l, permutation = zip(*list_vals)
- permutation = np.array(permutation)
- Uc = Uc[:, permutation]
- # And also rearrange the unitary and values so that they are decreasingly ordered
- return np.array(sorted_l), Uc
-
- v, l, ws = np.linalg.svd(N)
- w = np.transpose(np.conjugate(ws))
- rl = np.round(l, rounding)
-
- # Generate list with degenerancies
- result = []
- for k, g in groupby(rl):
- result.append(list(g))
-
- # Generate lists containing the columns that correspond to degenerancies
- kk = 0
- for k in result:
- for ind, j in enumerate(k): # pylint: disable=unused-variable
- k[ind] = kk
- kk = kk + 1
-
- # Generate the lists with the degenerate column subspaces
- vas = []
- was = []
- for i in result:
- vas.append(v[:, i])
- was.append(w[:, i])
-
- # Generate the matrices qs of the degenerate subspaces
- qs = []
- for i in range(len(result)):
- qs.append(sqrtm(np.transpose(vas[i]) @ was[i]))
-
- # Construct the Takagi unitary
- qb = block_diag(*qs)
-
- U = v @ np.conj(qb)
- return rl, U
-
+from thewalrus.symplectic import sympmat
+from thewalrus.decompositions import takagi
def graph_embed_deprecated(A, max_mean_photon=1.0, make_traceless=False, rtol=1e-05, atol=1e-08):
r"""Embed a graph into a Gaussian state.
@@ -140,7 +62,7 @@ def graph_embed_deprecated(A, max_mean_photon=1.0, make_traceless=False, rtol=1e
if make_traceless:
A = A - np.trace(A) * np.identity(n) / n
- s, U = takagi(A, tol=atol)
+ s, U = takagi(A)
sc = np.sqrt(1.0 + 1.0 / max_mean_photon)
vals = -np.arctanh(s / (s[0] * sc))
return vals, U
@@ -184,7 +106,7 @@ def graph_embed(A, mean_photon_per_mode=1.0, make_traceless=False, rtol=1e-05, a
scale = adj_scaling(A, n * mean_photon_per_mode)
A = scale * A
- s, U = takagi(A, tol=atol)
+ s, U = takagi(A)
vals = -np.arctanh(s)
return vals, U
@@ -903,172 +825,6 @@ def rectangular_compact(U, rtol=1e-12, atol=1e-12):
return _absorb_zeta(phases_temp)
-def williamson(V, tol=1e-11):
- r"""Williamson decomposition of positive-definite (real) symmetric matrix.
-
- See :ref:`williamson`.
-
- Note that it is assumed that the symplectic form is
-
- .. math:: \Omega = \begin{bmatrix}0&I\\-I&0\end{bmatrix}
-
- where :math:`I` is the identity matrix and :math:`0` is the zero matrix.
-
- See https://math.stackexchange.com/questions/1171842/finding-the-symplectic-matrix-in-williamsons-theorem/2682630#2682630
-
- Args:
- V (array[float]): positive definite symmetric (real) matrix
- tol (float): the tolerance used when checking if the matrix is symmetric: :math:`|V-V^T| \leq` tol
-
- Returns:
- tuple[array,array]: ``(Db, S)`` where ``Db`` is a diagonal matrix
- and ``S`` is a symplectic matrix such that :math:`V = S^T Db S`
- """
- (n, m) = V.shape
-
- if n != m:
- raise ValueError("The input matrix is not square")
-
- diffn = np.linalg.norm(V - np.transpose(V))
-
- if diffn >= tol:
- raise ValueError("The input matrix is not symmetric")
-
- if n % 2 != 0:
- raise ValueError("The input matrix must have an even number of rows/columns")
-
- n = n // 2
- omega = sympmat(n)
- vals = np.linalg.eigvalsh(V)
-
- for val in vals:
- if val <= 0:
- raise ValueError("Input matrix is not positive definite")
-
- Mm12 = sqrtm(np.linalg.inv(V)).real
- r1 = Mm12 @ omega @ Mm12
- s1, K = schur(r1)
- X = np.array([[0, 1], [1, 0]])
- I = np.identity(2)
- seq = []
-
- # In what follows I construct a permutation matrix p so that the Schur matrix has
- # only positive elements above the diagonal
- # Also the Schur matrix uses the x_1,p_1, ..., x_n,p_n ordering thus I use rotmat to
- # go to the ordering x_1, ..., x_n, p_1, ... , p_n
-
- for i in range(n):
- if s1[2 * i, 2 * i + 1] > 0:
- seq.append(I)
- else:
- seq.append(X)
-
- p = block_diag(*seq)
- Kt = K @ p
- s1t = p @ s1 @ p
- dd = xpxp_to_xxpp(s1t)
- perm_indices = xpxp_to_xxpp(np.arange(2 * n))
- Ktt = Kt[:, perm_indices]
- Db = np.diag([1 / dd[i, i + n] for i in range(n)] + [1 / dd[i, i + n] for i in range(n)])
- S = Mm12 @ Ktt @ sqrtm(Db)
- return Db, np.linalg.inv(S).T
-
-
-def bloch_messiah(S, tol=1e-10, rounding=9):
- r"""Bloch-Messiah decomposition of a symplectic matrix.
-
- See :ref:`bloch_messiah`.
-
- Decomposes a symplectic matrix into two symplectic unitaries and squeezing transformation.
- It automatically sorts the squeezers so that they respect the canonical symplectic form.
-
- Note that it is assumed that the symplectic form is
-
- .. math:: \Omega = \begin{bmatrix}0&I\\-I&0\end{bmatrix}
-
- where :math:`I` is the identity matrix and :math:`0` is the zero matrix.
-
- As in the Takagi decomposition, the singular values of N are considered
- equal if they are equal after np.round(values, rounding).
-
- If S is a passive transformation, then return the S as the first passive
- transformation, and set the the squeezing and second unitary matrices to
- identity. This choice is not unique.
-
- For more info see:
- https://math.stackexchange.com/questions/1886038/finding-euler-decomposition-of-a-symplectic-matrix
-
- Args:
- S (array[float]): symplectic matrix
- tol (float): the tolerance used when checking if the matrix is symplectic:
- :math:`|S^T\Omega S-\Omega| \leq tol`
- rounding (int): the number of decimal places to use when rounding the singular values
-
- Returns:
- tuple[array]: Returns the tuple ``(ut1, st1, vt1)``. ``ut1`` and ``vt1`` are symplectic orthogonal,
- and ``st1`` is diagonal and of the form :math:`= \text{diag}(s1,\dots,s_n, 1/s_1,\dots,1/s_n)`
- such that :math:`S = ut1 st1 v1`
- """
- (n, m) = S.shape
-
- if n != m:
- raise ValueError("The input matrix is not square")
- if n % 2 != 0:
- raise ValueError("The input matrix must have an even number of rows/columns")
-
- n = n // 2
- omega = sympmat(n)
- if np.linalg.norm(np.transpose(S) @ omega @ S - omega) >= tol:
- raise ValueError("The input matrix is not symplectic")
-
- if np.linalg.norm(np.transpose(S) @ S - np.eye(2 * n)) >= tol:
-
- u, sigma = polar(S, side="left")
- ss, uss = takagi(sigma, tol=tol, rounding=rounding)
-
- # Apply a permutation matrix so that the squeezers appear in the order
- # s_1,...,s_n, 1/s_1,...1/s_n
- perm = np.array(list(range(0, n)) + list(reversed(range(n, 2 * n))))
-
- pmat = np.identity(2 * n)[perm, :]
- ut = uss @ pmat
-
- # Apply a second permutation matrix to permute s
- # (and their corresonding inverses) to get the canonical symplectic form
- qomega = np.transpose(ut) @ (omega) @ ut
- st = pmat @ np.diag(ss) @ pmat
-
- # Identifying degenerate subspaces
- result = []
- for _k, g in groupby(np.round(np.diag(st), rounding)[:n]):
- result.append(list(g))
-
- stop_is = list(np.cumsum([len(res) for res in result]))
- start_is = [0] + stop_is[:-1]
-
- # Rotation matrices (not permutations) based on svd.
- # See Appendix B2 of Serafini's book for more details.
- u_list, v_list = [], []
-
- for start_i, stop_i in zip(start_is, stop_is):
- x = qomega[start_i:stop_i, n + start_i : n + stop_i].real
- u_svd, _s_svd, v_svd = np.linalg.svd(x)
- u_list = u_list + [u_svd]
- v_list = v_list + [v_svd.T]
-
- pmat1 = block_diag(*(u_list + v_list))
-
- st1 = pmat1.T @ pmat @ np.diag(ss) @ pmat @ pmat1
- ut1 = uss @ pmat @ pmat1
- v1 = np.transpose(ut1) @ u
-
- else:
- ut1 = S
- st1 = np.eye(2 * n)
- v1 = np.eye(2 * n)
-
- return ut1.real, st1.real, v1.real
-
def sun_compact(U, rtol=1e-12, atol=1e-12):
r"""Recursive factorization of unitary transfomations.
diff --git a/strawberryfields/ops.py b/strawberryfields/ops.py
index 68ff22505..fe05b47d3 100644
--- a/strawberryfields/ops.py
+++ b/strawberryfields/ops.py
@@ -26,6 +26,7 @@
import scipy.special as ssp
from thewalrus.symplectic import xxpp_to_xpxp
+from thewalrus.decompositions import blochmessiah, williamson
import strawberryfields as sf
import strawberryfields.program_utils as pu
@@ -2933,7 +2934,7 @@ class GaussianTransform(Decomposition):
transformation (:math:`S = O_1 Z`).
"""
- def __init__(self, S, vacuum=False, tol=1e-10):
+ def __init__(self, S, vacuum=False):
super().__init__([S])
self.ns = S.shape[0] // 2
self.vacuum = (
@@ -2954,7 +2955,7 @@ def __init__(self, S, vacuum=False, tol=1e-10):
self.U1 = X1 + 1j * P1
else:
# transformation is active, do Bloch-Messiah
- O1, smat, O2 = dec.bloch_messiah(S, tol=tol)
+ O1, smat, O2 = blochmessiah(S)
X1 = O1[:N, :N]
P1 = O1[N:, :N]
X2 = O2[:N, :N]
@@ -3072,7 +3073,7 @@ def __init__(self, V, r=None, decomp=True, tol=1e-6):
self.p_disp = r[self.ns :]
# needed only if decomposed
- th, self.S = dec.williamson(V, tol=tol)
+ th, self.S = williamson(V, rtol=tol)
self.pure = np.abs(np.linalg.det(V) - 1.0) < tol
self.nbar = 0.5 * (np.diag(th)[: self.ns] - 1.0)
diff --git a/strawberryfields/utils/random_numbers_matrices.py b/strawberryfields/utils/random_numbers_matrices.py
index 5062f2780..6b37f07af 100644
--- a/strawberryfields/utils/random_numbers_matrices.py
+++ b/strawberryfields/utils/random_numbers_matrices.py
@@ -110,12 +110,10 @@ def random_interferometer(N, real=False):
Returns:
array: random :math:`N\times N` unitary distributed with the Haar measure
"""
+ if N == 1:
+ if real:
+ return np.array([[2 * (np.random.binomial(1, 0.5) - 0.5)]])
+ return np.array([[np.exp(1j * 2 * np.pi * np.random.rand())]])
if real:
- z = np.random.randn(N, N)
- else:
- z = randnc(N, N) / np.sqrt(2.0)
- q, r = sp.linalg.qr(z)
- d = np.diagonal(r)
- ph = d / np.abs(d)
- U = np.multiply(q, ph, q)
- return U
+ return sp.stats.ortho_group.rvs(N)
+ return sp.stats.unitary_group.rvs(N)
diff --git a/tests/frontend/test_decompositions.py b/tests/frontend/test_decompositions.py
index 9df20787d..0a4cb341a 100644
--- a/tests/frontend/test_decompositions.py
+++ b/tests/frontend/test_decompositions.py
@@ -22,6 +22,7 @@
import numpy as np
import scipy as sp
from scipy.linalg import qr, block_diag
+from thewalrus.symplectic import sympmat as omega
from strawberryfields import decompositions as dec
from strawberryfields.utils import random_interferometer as haar_measure
@@ -29,63 +30,6 @@
N_SAMPLES = 10
-def omega(n):
- """Returns the symplectic matrix for n modes"""
- idm = np.identity(n)
- O = np.concatenate(
- (
- np.concatenate((0 * idm, idm), axis=1),
- np.concatenate((-idm, 0 * idm), axis=1),
- ),
- axis=0,
- )
- return O
-
-
-class TestTakagi:
- """Takagi decomposition tests"""
-
- def test_square_validation(self):
- """Test that the takagi decomposition raises exception if not square"""
- A = np.random.random([4, 5]) + 1j * np.random.random([4, 5])
- with pytest.raises(ValueError, match="matrix must be square"):
- dec.takagi(A)
-
- def test_symmetric_validation(self):
- """Test that the takagi decomposition raises exception if not symmetric"""
- A = np.random.random([5, 5]) + 1j * np.random.random([5, 5])
- with pytest.raises(ValueError, match="matrix is not symmetric"):
- dec.takagi(A)
-
- def test_random_symm(self, tol):
- """Verify that the Takagi decomposition, applied to a random symmetric
- matrix, produced a decomposition that can be used to reconstruct the matrix."""
- A = np.random.random([6, 6]) + 1j * np.random.random([6, 6])
- A += A.T
- rl, U = dec.takagi(A)
- res = U @ np.diag(rl) @ U.T
- assert np.allclose(res, A, atol=tol, rtol=0)
-
- def test_real_degenerate(self):
- """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."""
- g = nx.balanced_tree(2, 4)
- a = nx.to_numpy_array(g)
- rl, U = dec.takagi(a)
- assert np.allclose(U @ U.conj().T, np.eye(len(a)))
- assert np.allclose(U @ np.diag(rl) @ U.T, a)
-
- def test_zeros(self):
- """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 = dec.takagi(a)
- assert np.allclose(rl, np.zeros(dim))
- assert np.allclose(U, np.eye(dim))
-
-
class TestGraphEmbed:
"""graph_embed tests"""
@@ -546,258 +490,6 @@ def test_decomposition(self, U, tol):
assert np.allclose(U, Uout, atol=tol, rtol=0)
-class TestWilliamsonDecomposition:
- """Tests for the Williamson decomposition"""
-
- @pytest.fixture
- def create_cov(self, hbar, tol):
- """create a covariance state for use in testing.
-
- Args:
- nbar (array[float]): vector containing thermal state values
-
- Returns:
- tuple: covariance matrix and symplectic transform
- """
-
- def _create_cov(nbar):
- """wrapped function"""
- n = len(nbar)
- O = omega(n)
-
- # initial vacuum state
- cov = np.diag(2 * np.tile(nbar, 2) + 1) * hbar / 2
-
- # interferometer 1
- U1 = haar_measure(n)
- S1 = np.vstack([np.hstack([U1.real, -U1.imag]), np.hstack([U1.imag, U1.real])])
-
- # squeezing
- r = np.log(0.2 * np.arange(n) + 2)
- Sq = block_diag(np.diag(np.exp(-r)), np.diag(np.exp(r)))
-
- # interferometer 2
- U2 = haar_measure(n)
- S2 = np.vstack([np.hstack([U2.real, -U2.imag]), np.hstack([U2.imag, U2.real])])
-
- # final symplectic
- S_final = S2 @ Sq @ S1
-
- # final covariance matrix
- cov_final = S_final @ cov @ S_final.T
-
- # check valid symplectic transform
- assert np.allclose(S_final.T @ O @ S_final, O)
-
- # check valid state
- eigs = np.linalg.eigvalsh(cov_final + 1j * (hbar / 2) * O)
- eigs[np.abs(eigs) < tol] = 0
- assert np.all(eigs >= 0)
-
- if np.allclose(nbar, 0):
- # check pure
- assert np.allclose(np.linalg.det(cov_final), (hbar / 2) ** (2 * n))
- else:
- # check not pure
- assert not np.allclose(np.linalg.det(cov_final), (hbar / 2) ** (2 * n))
-
- return cov_final, S_final
-
- return _create_cov
-
- def test_square_validation(self):
- """Test that the graph_embed decomposition raises exception if not square"""
- A = np.random.random([4, 5]) + 1j * np.random.random([4, 5])
- with pytest.raises(ValueError, match="matrix is not square"):
- dec.williamson(A)
-
- def test_symmetric_validation(self):
- """Test that the graph_embed decomposition raises exception if not symmetric"""
- A = np.random.random([5, 5]) + 1j * np.random.random([5, 5])
- with pytest.raises(ValueError, match="matrix is not symmetric"):
- dec.williamson(A)
-
- def test_even_validation(self):
- """Test that the graph_embed decomposition raises exception if not even number of rows"""
- A = np.random.random([5, 5]) + 1j * np.random.random([5, 5])
- A += A.T
- with pytest.raises(ValueError, match="must have an even number of rows/columns"):
- dec.williamson(A)
-
- def test_positive_definite_validation(self):
- """Test that the graph_embed decomposition raises exception if not positive definite"""
- A = np.diag([-2, 0.1, 2, 3])
- with pytest.raises(ValueError, match="matrix is not positive definite"):
- dec.williamson(A)
-
- def test_vacuum_state(self, tol, hbar):
- """Test vacuum state"""
- V = np.identity(4)
- Db, S = dec.williamson(V)
- assert np.allclose(Db, np.identity(4), atol=tol, rtol=0)
- assert np.allclose(S, np.identity(4), atol=tol, rtol=0)
-
- def test_pure_state(self, create_cov, hbar, tol):
- """Test pure state"""
- n = 3
- O = omega(n)
-
- cov, _ = create_cov(np.zeros([n]))
-
- Db, S = dec.williamson(cov)
- nbar = np.diag(Db) / hbar - 0.5
-
- # check decomposition is correct
- assert np.allclose(S @ Db @ S.T, cov, atol=tol, rtol=0)
- # check nbar = 0
- assert np.allclose(nbar, 0, atol=tol, rtol=0)
- # check S is symplectic
- assert np.allclose(S @ O @ S.T, O, atol=tol, rtol=0)
-
- def test_mixed_state(self, create_cov, hbar, tol):
- """Test mixed state"""
- n = 3
- O = omega(n)
- nbar_in = np.abs(np.random.random(n))
-
- cov, _ = create_cov(nbar_in)
-
- Db, S = dec.williamson(cov)
- nbar = np.diag(Db) / hbar - 0.5
-
- # check decomposition is correct
- assert np.allclose(S @ Db @ S.T, cov, atol=tol, rtol=0)
- # check nbar
- assert np.allclose(sorted(nbar[:n]), sorted(nbar_in), atol=tol, rtol=0)
- # check S is symplectic
- assert np.allclose(S @ O @ S.T, O, atol=tol, rtol=0)
-
-
-class TestBlochMessiahDecomposition:
- """Tests for the Bloch-Messiah decomposition"""
-
- @pytest.fixture
- def create_transform(self):
- """create a symplectic transform for use in testing.
-
- Args:
- n (int): number of modes
- passive (bool): whether transform should be passive or not
-
- Returns:
- array: symplectic matrix
- """
-
- def _create_transform(n, passive=True):
- """wrapped function"""
- O = omega(n)
-
- # interferometer 1
- U1 = haar_measure(n)
- S1 = np.vstack([np.hstack([U1.real, -U1.imag]), np.hstack([U1.imag, U1.real])])
-
- Sq = np.identity(2 * n)
- if not passive:
- # squeezing
- r = np.log(0.2 * np.arange(n) + 2)
- Sq = block_diag(np.diag(np.exp(-r)), np.diag(np.exp(r)))
-
- # interferometer 2
- U2 = haar_measure(n)
- S2 = np.vstack([np.hstack([U2.real, -U2.imag]), np.hstack([U2.imag, U2.real])])
-
- # final symplectic
- S_final = S2 @ Sq @ S1
-
- # check valid symplectic transform
- assert np.allclose(S_final.T @ O @ S_final, O)
- return S_final
-
- return _create_transform
-
- def test_square_validation(self):
- """Test raises exception if not square"""
- A = np.random.random([4, 5]) + 1j * np.random.random([4, 5])
- with pytest.raises(ValueError, match="matrix is not square"):
- dec.bloch_messiah(A)
-
- def test_symmplectic(self):
- """Test raises exception if not symmetric"""
- A = np.random.random([6, 6]) + 1j * np.random.random([6, 6])
- A += A.T
- with pytest.raises(ValueError, match="matrix is not symplectic"):
- dec.bloch_messiah(A)
-
- def test_even_validation(self):
- """Test raises exception if not even number of rows"""
- A = np.random.random([5, 5]) + 1j * np.random.random([5, 5])
- A += A.T
- with pytest.raises(ValueError, match="must have an even number of rows/columns"):
- dec.bloch_messiah(A)
-
- def test_identity(self, tol):
- """Test identity"""
- n = 2
- S_in = np.identity(2 * n)
- O1, S, O2 = dec.bloch_messiah(S_in)
-
- assert np.allclose(O1 @ O2, np.identity(2 * n), atol=tol, rtol=0)
- assert np.allclose(S, np.identity(2 * n), atol=tol, rtol=0)
-
- # test orthogonality
- assert np.allclose(O1.T, O1, atol=tol, rtol=0)
- assert np.allclose(O2.T, O2, atol=tol, rtol=0)
-
- # test symplectic
- O = omega(n)
- assert np.allclose(O1 @ O @ O1.T, O, atol=tol, rtol=0)
- assert np.allclose(O2 @ O @ O2.T, O, atol=tol, rtol=0)
-
- def test_passive_transform(self, create_transform, tol):
- """Test passive transform has no squeezing.
- Note: this test also tests the case with degenerate symplectic values"""
- n = 3
- S_in = create_transform(3, passive=True)
- O1, S, O2 = dec.bloch_messiah(S_in)
-
- # test decomposition
- assert np.allclose(O1 @ S @ O2, S_in, atol=tol, rtol=0)
-
- # test no squeezing
- assert np.allclose(O1 @ O2, S_in, atol=tol, rtol=0)
- assert np.allclose(S, np.identity(2 * n), atol=tol, rtol=0)
-
- # test orthogonality
- assert np.allclose(O1.T @ O1, np.identity(2 * n), atol=tol, rtol=0)
- assert np.allclose(O2.T @ O2, np.identity(2 * n), atol=tol, rtol=0)
-
- # test symplectic
- O = omega(n)
- # TODO: BUG:
- # 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)
-
- def test_active_transform(self, create_transform, tol):
- """Test passive transform with squeezing"""
- n = 3
- S_in = create_transform(3, passive=False)
- O1, S, O2 = dec.bloch_messiah(S_in)
-
- # test decomposition
- assert np.allclose(O1 @ S @ O2, S_in, atol=tol, rtol=0)
-
- # test orthogonality
- assert np.allclose(O1.T @ O1, np.identity(2 * n), atol=tol, rtol=0)
- assert np.allclose(O2.T @ O2, np.identity(2 * n), atol=tol, rtol=0)
-
- # test symplectic
- O = omega(n)
- 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)
-
-
class TestSUnFactorization:
"""tests for the SU(n) factorization"""
diff --git a/tests/frontend/test_ops_decompositions.py b/tests/frontend/test_ops_decompositions.py
index 773832cf7..bedef1958 100644
--- a/tests/frontend/test_ops_decompositions.py
+++ b/tests/frontend/test_ops_decompositions.py
@@ -19,7 +19,8 @@
import numpy as np
from thewalrus.quantum import Amat
-
+from thewalrus.symplectic import expand, rotation, beam_splitter, squeezing, two_mode_squeezing
+from thewalrus.decompositions import blochmessiah
import strawberryfields as sf
from strawberryfields.parameters import par_evaluate, FreeParameter
from strawberryfields import decompositions as dec
@@ -81,11 +82,8 @@ def _rotation(phi, mode, num_modes):
Returns:
array[float]: transformation matrix
"""
- c = np.cos(phi)
- s = np.sin(phi)
- S = np.array([[c, -s], [s, c]])
- return expand(S, mode, num_modes)
+ return expand(rotation(phi), mode, num_modes)
def _squeezing(r, phi, mode, num_modes):
@@ -100,14 +98,8 @@ def _squeezing(r, phi, mode, num_modes):
Returns:
array: symplectic transformation matrix
"""
- cp = np.cos(phi)
- sp = np.sin(phi)
- ch = np.cosh(r)
- sh = np.sinh(r)
-
- S = np.array([[ch - cp * sh, -sp * sh], [-sp * sh, ch + cp * sh]])
- return expand(S, mode, num_modes)
+ return expand(squeezing(r, phi), mode, num_modes)
def _two_mode_squeezing(r, phi, modes, num_modes):
@@ -122,21 +114,8 @@ def _two_mode_squeezing(r, phi, modes, num_modes):
Returns:
array: symplectic transformation matrix
"""
- cp = np.cos(phi)
- sp = np.sin(phi)
- ch = np.cosh(r)
- sh = np.sinh(r)
-
- S = np.array(
- [
- [ch, cp * sh, 0, sp * sh],
- [cp * sh, ch, sp * sh, 0],
- [0, sp * sh, ch, -cp * sh],
- [sp * sh, 0, -cp * sh, ch],
- ]
- )
- return expand(S, modes, num_modes)
+ return expand(two_mode_squeezing(r, phi), modes, num_modes)
def _beamsplitter(theta, phi, modes, num_modes):
@@ -151,21 +130,8 @@ def _beamsplitter(theta, phi, modes, num_modes):
Returns:
array[float]: transformation matrix
"""
- cp = np.cos(phi)
- sp = np.sin(phi)
- ct = np.cos(theta)
- st = np.sin(theta)
-
- S = np.array(
- [
- [ct, -cp * st, 0, -st * sp],
- [cp * st, ct, -st * sp, 0],
- [0, st * sp, ct, -cp * st],
- [st * sp, 0, cp * st, ct],
- ]
- )
- return expand(S, modes, num_modes)
+ return expand(beam_splitter(theta, phi), modes, num_modes)
class TestDecompositions:
@@ -513,7 +479,7 @@ def test_decomposition_active(self, tol):
n = 3
S = random_symplectic(n, passive=False)
- O1, Sq, O2 = dec.bloch_messiah(S)
+ O1, Sq, O2 = blochmessiah(S)
X1 = O1[:n, :n]
P1 = O1[n:, :n]
X2 = O2[:n, :n]
@@ -594,7 +560,7 @@ def test_active_on_vacuum(self, tol):
n = 3
S = random_symplectic(n, passive=False)
- O1, _, _ = dec.bloch_messiah(S)
+ O1, _, _ = blochmessiah(S)
X1 = O1[:n, :n]
P1 = O1[n:, :n]
# X2 = O2[:n, :n]
diff --git a/tests/frontend/test_utils.py b/tests/frontend/test_utils.py
index 05562f02c..06c0dd575 100644
--- a/tests/frontend/test_utils.py
+++ b/tests/frontend/test_utils.py
@@ -267,14 +267,10 @@ def test_odd_cat_state(self, a, cutoff, tol):
# ===================================================================================
+@pytest.mark.parametrize("modes", [1, 2, 3])
class TestRandomMatrices:
"""Unit tests for random matrices"""
- @pytest.fixture
- def modes(self):
- """Number of modes to use when creating matrices"""
- return 3
-
@pytest.mark.parametrize("pure_state", [True, False])
@pytest.mark.parametrize("block_diag", [True, False])
def test_random_covariance_square(self, modes, hbar, pure_state, block_diag):