Skip to content

inverse_problem

qoptcraft.evolution.inverse_problem

scattering_from_unitary(unitary, modes, photons)

Retrieve the linear optical scattering matrix from a unitary using Theorem 2 of PRA 100, 022301 (2019).

Parameters:

Name Type Description Default
unitary NDArray

unitary matrix.

 required
modes int

number of modes in the optical system.

 required
photons int

number of photons.

 required

Returns:

Name Type Description
NDArray NDArray

optical scattering matrix that maps to the given unitary.
Source code in qoptcraft/evolution/inverse_problem.py 
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
68
69
70
71
def scattering_from_unitary(unitary: NDArray, modes: int, photons: int) -> NDArray:
    """Retrieve the linear optical scattering matrix from a unitary using Theorem 2
    of PRA 100, 022301 (2019).

    Args:
        unitary (NDArray): unitary matrix.
        modes (int): number of modes in the optical system.
        photons (int): number of photons.

    Returns:
        NDArray: optical scattering matrix that maps to the given unitary.
    """
    coefs = _solution_coefs(unitary, modes, photons)
    basis, _ = get_algebra_basis(modes, photons)

    basis_array = np.array(basis)

    def adjoint(basis_matrix):
        for i in range(len(basis)):
            if (basis_matrix == basis[i]).all():
                index = i
                break
        return np.einsum("k,kij->ij", coefs[index, :], basis_array)

    identity = np.eye(modes)

    def get_nonzero_element():
        for j in range(modes):
            basis_matrix = sym_matrix(j, j, modes)
            for l in range(modes):
                ad_matrix = adjoint(basis_matrix)
                exp_val = np.einsum("i,ij,j->", identity[l, :], ad_matrix, identity[l, :])
                if not np.isclose(exp_val, 0, rtol=1e-5, atol=1e-8):
                    j0, l0 = j, l
                    coef = np.sqrt(-1j * exp_val)
                    return j0, l0, coef
        raise ValueError("Nonzero element not found.")

    j0, l0, coef = get_nonzero_element()
    scattering = np.zeros((modes, modes), dtype=np.complex128)
    for l in range(modes):
        for j in range(modes):
            ad_sym = adjoint(sym_matrix(j, j0, modes))
            sym_term = np.einsum("i,ij,j->", identity[l, :], ad_sym, identity[l0, :])
            if j != j0:
                ad_antisym = adjoint(antisym_matrix(j, j0, modes))
                antisym_term = np.einsum(
                    "i,ij,j->",
                    identity[l, :],
                    ad_antisym,
                    identity[l0, :],
                )
            else:
                antisym_term = 0
            scattering[l, j] = (antisym_term - 1j * sym_term) / coef
    return scattering