Banded unitaries

.github/CHANGELOG.md

@@ -2,6 +2,8 @@
 
 ### New features
 
+* Adds the function `random_banded_interferometer` to generate unitary matrices with a given bandwidth. [#208](https://github.com/XanaduAI/thewalrus/pull/208)
+
 ### Improvements
 
 ### Bug fixes
@@ -12,6 +14,8 @@
 
 This release contains contributions from (in alphabetical order):
 
+Nicolas Quesada
+
 ---
 
 # Version 0.14.0

thewalrus/random.py

@@ -117,3 +117,53 @@ def random_interferometer(N, real=False):
     ph = d / np.abs(d)
     U = np.multiply(q, ph, q)
     return U
+
+
+def random_block_interferometer(N, top_one=True, real=False):
+    r"""Random interferometer with blocks of at most size 2.
+
+    Args:
+        N (int): number of modes
+        top_one (bool): if true places a `1\times1` interferometer in the top-left most block
+        real (bool): return a random real orthogonal matrix
+
+    Returns:
+        array: random :math:`N\times N` unitary with the specified block structure
+    """
+    if N % 2 == 0:
+        if top_one is True:
+            u2s = [random_interferometer(2, real=real) for i in range(((N // 2) - 1))]
+            u0 = random_interferometer(1, real=real)
+            u1 = random_interferometer(1, real=real)
+            return sp.linalg.block_diag(u0, *u2s, u1)
+
+        u2s = [random_interferometer(2, real=real) for i in range(N // 2)]
+        return sp.linalg.block_diag(*u2s)
+
+    u2s = [random_interferometer(2, real=real) for i in range((N - 1) // 2)]
+    u0 = random_interferometer(1, real=real)
+    if top_one is True:
+        return sp.linalg.block_diag(u0, *u2s)
+    return sp.linalg.block_diag(*u2s, u0)
+
+def random_banded_interferometer(N, w, top_one_init=True, real=False):
+    r"""Generates a banded unitary matrix
+
+    Args:
+        N (int): number of modes
+        w (int): bandwidth
+        top_one_init (bool): if true places a `1\times1` interferometer in the top-left-most block of the first matrix in the product
+        real (bool): return a random real orthogonal matrix
+
+    Returns:
+        array: random :math:`N\times N` unitary with the specified block structure
+    """
+    if N < w + 1:
+        raise ValueError("The bandwidth can be at most one minus the size of the matrix.")
+    if N == w + 1:
+        return random_interferometer(N, real=real)
+    U = sp.linalg.block_diag(*[random_interferometer(1, real=real) for _ in range(N)])
+    for _ in range(w):
+        U = U @ random_block_interferometer(N, top_one=top_one_init, real=real)
+        top_one_init = not top_one_init
+    return U

thewalrus/tests/test_random.py

@@ -0,0 +1,65 @@
+# Copyright 2019 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.
+# You may obtain a copy of the License at
+
+#     http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+"""Random tests"""
+import pytest
+
+import numpy as np
+from thewalrus.random import random_block_interferometer, random_banded_interferometer
+
+
+def bandwidth(A):
+    """Calculates the upper bandwidth of the matrix A
+
+    Args:
+        A (array): input matrix
+
+    Returns:
+        (int): bandwidth of matrix
+    """
+    n, _ = A.shape
+    for i in range(n):
+        vali = np.diag(A, i)
+        if np.allclose(vali, 0):
+            return i - 1
+    return n - 1
+
+
+@pytest.mark.parametrize("n", [5, 7, 8, 9])
+@pytest.mark.parametrize("top_one", [True, False])
+@pytest.mark.parametrize("real", [True, False])
+def test_random_block(n, top_one, real):
+    """Test that random_block_interferometer produces a unitary with the right structure"""
+    U = random_block_interferometer(n, top_one=top_one, real=real)
+    assert np.allclose(U @ U.T.conj(), np.identity(n))
+    if top_one:
+        assert np.allclose(U[0, 1], 0)
+        assert np.allclose(U[1, 0], 0)
+
+
+@pytest.mark.parametrize("n", [5, 7, 8, 9])
+@pytest.mark.parametrize("top_one_init", [True, False])
+@pytest.mark.parametrize("real", [True, False])
+def test_random_banded(n, top_one_init, real):
+    """Test that random_banded_interferometer produces a unitary with the right structure"""
+    for w in range(n):
+        U = random_banded_interferometer(n, w, top_one_init=top_one_init, real=real)
+        assert np.allclose(U @ U.T.conj(), np.identity(n))
+        assert w == bandwidth(U)
+
+def test_wrong_bandwidth():
+    """Test that the correct error is raised if w > n-1"""
+    n = 10
+    w = 10
+    with pytest.raises(ValueError, match="The bandwidth can be at most one minus the size of the matrix."):
+        random_banded_interferometer(n, w)