Skip to content

clemens_decomp

qoptcraft.optical_elements.clemens_decomp

Decompose a unitary into beamsplitters following Clemens et al.

clemens_decomposition(unitary)

Given a unitary matrix calculates the Clemens et al. decompositon into beam splitters and phase shifters:

D = L_n ... L_1 · U · R_1.inv ... R_n.inv => => U = L_1.inv ... L_n.inv · D · R_n ... R_1

Parameters:

Name Type Description Default
unitary NDArray

unitary matrix to decompose.

 required

Returns:

Type Description
tuple[list[NDArray], NDArray, list[NDArray]]

list[NDArray]: list of matrices that decompose U in the order of the decomposition
References

[1] Clements et al., "An Optimal Design for Universal Multiport Interferometers" arXiv, 2007. https://arxiv.org/pdf/1603.08788.pdf

Source code in qoptcraft/optical_elements/clemens_decomp.py 
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
def clemens_decomposition(unitary: NDArray) -> tuple[list[NDArray], NDArray, list[NDArray]]:
    """Given a unitary matrix calculates the Clemens et al. decompositon
    into beam splitters and phase shifters:

    D = L_n ... L_1 · U · R_1.inv ... R_n.inv  =>
    => U = L_1.inv ... L_n.inv · D · R_n ... R_1

    Args:
        unitary (NDArray): unitary matrix to decompose.

    Returns:
        list[NDArray]: list of matrices that decompose U in the order of the decomposition

    References:
        [1] Clements et al., "An Optimal Design for Universal Multiport
            Interferometers" arXiv, 2007. https://arxiv.org/pdf/1603.08788.pdf
    """

    assert unitary.shape[0] == unitary.shape[1], "The matrix is not square"
    dim = unitary.shape[0]

    right_list = []
    left_list = []

    for i in range(1, dim):
        if (i % 2) == 1:
            for j in range(i):
                # substract 1 to all indices to match the paper with python arrays
                row = dim - j - 1
                col = i - j - 1
                mode_1 = i - j - 1  # mode_1 must equal col
                mode_2 = i - j

                # We calculate the beamsplitter angles phi y theta
                angle, shift = _solve_angles(unitary, mode_1, mode_2, row, col, is_odd=True)
                R = beam_splitter(angle, shift, dim, mode_1, mode_2, convention="clemens")
                unitary = unitary @ R.conj().T
                right_list.append(R)
        else:
            for j in range(1, i + 1):
                # substract 1 to all indices to match the paper with python arrays
                row = dim + j - i - 1
                col = j - 1
                mode_1 = dim + j - i - 2
                mode_2 = dim + j - i - 1  # mode_2 must equal row

                # We calculate the beamsplitter angles phi y theta
                angle, shift = _solve_angles(unitary, mode_1, mode_2, row, col, is_odd=False)
                left = beam_splitter(angle, shift, dim, mode_1, mode_2, convention="clemens")
                unitary = left @ unitary
                left_list.append(left.conj().T)
    diag = unitary

    right_list.reverse()  # save as [R_n, ..., R_1]
    # left_list = [L_1.inv, ... L_n.inv]

    return left_list, diag, right_list