Skip to content

reck_decomp

qoptcraft.optical_elements.reck_decomp

Decompose a unitary into beamsplitters following Reck et al.

reck_decomposition(unitary)

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

D = U · BS_1... BS_n => U = D · BS_n.inv ... BS_1.inv

Parameters:

Name Type Description Default
unitary NDArray

unitary matrix to decompose.

 required

Returns:

Type Description
tuple[NDArray, list[NDArray]]

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

[1] Reck et al., "Experimental realization of any discrete unitary operator" Phys. Rev. Lett. 73, 58, 1994. https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.73.58

Source code in qoptcraft/optical_elements/reck_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
def reck_decomposition(unitary: NDArray) -> tuple[NDArray, list[NDArray]]:
    """Given a unitary matrix calculates the Reck et al. decompositon
    into beam splitters and phase shifters:

    D = U · BS_1... BS_n  =>  U = D · BS_n.inv ... BS_1.inv

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

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

    References:
        [1] Reck et al., "Experimental realization of any discrete unitary operator"
            Phys. Rev. Lett. 73, 58, 1994. https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.73.58
    """

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

    bs_list = []

    for row in range(dim - 1, 0, -1):
        for col in range(row - 1, -1, -1):
            # Eliminate (row, col) element of U with a beamsplitter T_(row, col)

            # We calculate the beamsplitter angles phi y theta
            angle, shift = _solve_angles(unitary, row, col)
            right = beam_splitter(angle, shift, dim, mode_1=row, mode_2=col, convention="reck")
            unitary = unitary @ right
            bs_list.append(right.conj().T)
    diag = unitary
    bs_list.reverse()
    return diag, bs_list