Quantum pattern recognition in photonic circuits

Quantum information processing deals with the manipulation of quantum states in order to perform quantum informational tasks such as quantum algorithms [1], quantum error correction [2], quantum cryptography [3], or quantum teleportation [4]. It is known that quantum information has the potential to outperform classical information protocols [5–9]. However, since a quantum system is modified after measurement, extracting arbitrary information from a quantum state may require many copies. In general, quantum tomography (QT) can be employed for the full reconstruction of a quantum state or a quantum operator. In a nutshell, QT implies the measurement of the expectation values of several operators or the use of a mutually unbiased basis [10–12], which is in general a hard experimental task. In the last years, some efficient protocols have been proposed assisted by machine learning (ML) algorithms, although it may still be tricky depending on the dimension of the quantum system [13–18].

Quantum information can be encoded, decoded, and manipulated in a variety of physical systems, for example, photonic systems [19], solid state devices [20], trapped ions [21, 22], and superconducting architectures [23, 24]. Photonic platforms present several advantages for quantum information protocols, such as long coherence times and full connectivity, allowing long-distance quantum communication and quantum key distribution, among many other achievements [25–27]. Different photonic degrees of freedom, including polarization, spectral, spatial, and temporal modes can be used to encode information [19], providing a huge variety of experimental resources for many quantum information tasks.

One of the most intriguing experimental resources for photonic platform is boson sampling, which is a model of non-universal quantum computation. It consists in measuring or sampling the quantity of photons in a distribution produced by a linear interference device given an initial photonic Fock state [28]. It is a complex problem that cannot be efficiently simulated on conventional computers. However, with accessible experimental requirements, it can be simulated on linear optical quantum computers. It is attractive that the experimental setup of boson sampling only requires single-photon sources, photodetectors, and linear optical elements, i.e. beam splitters and phase shifters [29]. Such a feasibility has encouraged and inspired many research teams for lab implementations [30–34]. Particularly, recent experiments have shown the potential of photonic technologies to obtain quantum advantage for scientific relevant problems; for instance, quantum supremacy of Gaussian boson sampling with 76 photons [7], and the so-called timestamp boson sampling, that brings to light the potential of circuits to get memory effects, thus paving the way for the implementation of neuromorphic quantum computing [35, 36].

On a separate note, it is known that ML algorithms help to extract information from large and complex data sets, encompassing different techniques with sound mathematical grounds [37, 38]. The information is extracted by means of learning algorithms, that teach models which usually have no a priori knowledge of the problem. The increasing availability and size of data sets has spread the use of ML [39, 40], oftentimes involving reliable applications at academic and commercial levels.

The main goal of this manuscript is to show a simple and experimentally feasible photonic architecture capable to produce photonic patterns which assisted by ML algorithms could give an estimation about quantum features hard to calculate experimentally. To this end, we consider whether the full tomography of a photonic mode and two-mode entanglement can be extracted via sampling the output distribution of photons. In order to obtain a relation between the distribution probability of the output state and the information of the unknown state, we propose a particular linear optical circuit with four spatial modes, and calculate the corresponding permanent of the submatrix of the unitary matrix. We make use of a data set of probability patterns, that corresponds with the photon number distributions of output states. Then, we change the parameters of the unknown state to build a supervised learning algorithm for estimating the state of a new probability pattern, distinct from the training one.

The rest of the article is organized as follows. In section 2, we briefly review the most important aspects of boson sampling. Sections 3 and 4 introduce our four-mode optical circuit for state estimation of a photonic-mode, and entanglement estimation for a bipartite system, respectively. The results of the ML approach for state estimation and entanglement estimation are shown in section 5, ending up the paper with the conclusions of the work in section 6.

A boson-sampling experiment consists in measuring the photon number probability produced by the interference of N photonic states, usually indistinguishable single-photon states, via an M-mode linear network. The distribution can be obtained by computing permanents of the submatrix derived from the unitary transformation matrix of the network [41]. The calculation of the permanent for an n × n complex matrix is a #P-hard task for classical computers and #P-complete for a (0, 1)-matrix [42], which means that the simulation of a boson-sampling experiment by classical devices is inefficient.

Without loss of generality, the output state of an M-mode optical circuit after the interference of an N-photon state can be calculated as follows. First, the input state is given by

$$\begin{equation}\vert {\psi }_{\text{in}}\rangle =\left\vert {I}_{1},{I}_{2},\dots ,{I}_{M}\right\rangle =\left(\prod\limits _{k}\frac{{\hat{a}}_{k}^{{\dagger}{I}_{k}}}{\sqrt{{I}_{k}!}}\right)\left\vert 0\right\rangle ,\end{equation} \tag{ 1 }$$

where ${\hat{a}}_{k}^{{\dagger}}$ is the photon creation operator for the kth mode, I_k is the number of photons in the kth mode, and ∑_k I_k = N. The input and output states are related by a unitary transformation,

$$\begin{equation}\left\vert {\psi }_{\text{out}}\right\rangle =\mathbb{U}\left\vert {\psi }_{\text{in}}\right\rangle =\mathbb{U}\left(\prod\limits _{k}\frac{{\hat{a}}_{k}^{{\dagger}{I}_{k}}}{\sqrt{{I}_{k}!}}\right){\mathbb{U}}^{{\dagger}}\mathbb{U}\left\vert 0\right\rangle ,\end{equation} \tag{ 2 }$$

where $\mathbb{U}\left\vert 0\right\rangle =\left\vert 0\right\rangle$ because a linear optical circuit preserves the photon number. Now, the operator transformation reads

$$\begin{equation}\mathbb{U}{\hat{a}}_{k}^{{\dagger}}{\mathbb{U}}^{{\dagger}}=\sum\limits _{j}\hspace{2pt}{u}_{k,j}\enspace {a}_{j}^{{\dagger}}.\end{equation} \tag{ 3 }$$

Again, as a linear optical circuit preserves the number of photons, then u_k,j defines a unitary matrix U, which represents a superoperator that acts over the space of the creation operators. Then, the output state reads

$$\begin{equation}\left\vert {\psi }_{\text{out}}\right\rangle =\left[\prod\limits _{k}\frac{{(U{\hat{a}}_{k}^{{\dagger}})}^{{I}_{k}}}{\sqrt{{I}_{k}!}}\right]\left\vert 0\right\rangle =\sum\limits _{O}\hspace{2pt}{\gamma }_{O}\underset{k}{\bigotimes}\left\vert {O}_{k}\right\rangle =\sum\limits _{O}\hspace{2pt}{\gamma }_{O}\left\vert {\psi }_{\text{O}}\right\rangle ,\end{equation} \tag{ 4 }$$

where O is a photon-mode configuration containing N photons, γ_O is the superposition factor of the configuration O, and $\left\vert {O}_{k}\right\rangle$ is the photon number state for the kth-output mode in the Oth configuration. The probability of measuring configuration O is given by P_O = |γ_O |². The probability P_O is given by

$$\begin{equation}{P}_{O}=\frac{{\left\vert \text{Per}[{{\Lambda}}_{\vert {\psi }_{\text{I}}\rangle ,\vert {\psi }_{\text{O}}\rangle }]\right\vert }^{2}}{{I}_{1}!\dots {I}_{M}!{O}_{1}!\dots {O}_{M}!},\end{equation} \tag{ 5 }$$

where Per[⋅] is the permanent and ${{\Lambda}}_{\vert {\psi }_{\text{I}}\rangle ,\vert {\psi }_{\text{O}}\rangle }$ is an N × N matrix, which can be obtained from the elements u_jk of that defines U [29]. The matrix ${{\Lambda}}_{\vert {\psi }_{\text{I}}\rangle ,\vert {\psi }_{\text{O}}\rangle }$ reads

$$\begin{align}{({{\Lambda}}_{\vert {\psi }_{\text{I}}\rangle ,\vert {\psi }_{\text{O}}\rangle })}_{j,k}={U}_{p,q}={u}_{p,q}& \enspace \Longleftrightarrow\enspace {S}_{p-1}^{\text{O}}+1\leqslant j\leqslant {S}_{p}^{\text{O}},\\ & {S}_{q-1}^{\text{I}}+1\leqslant k\leqslant {S}_{q}^{\text{I}},\end{align} \tag{ 6 }$$

where

$$\begin{equation}{S}_{\ell }^{\text{O}}=\sum\limits _{i=1}^{\ell }{O}_{i}\quad \wedge \quad {S}_{\ell }^{\text{I}}=\sum\limits _{i=1}^{\ell }{I}_{i}\end{equation} \tag{ 7 }$$

are the number of output and input photons until mode ℓ, respectively, being ${S}_{0}^{\text{I}(\text{O})}=0$. It should be borne in mind that, for any linear network with M modes, we have that ${S}_{M}^{\text{I}}={S}_{M}^{\text{O}}=N$. For example, if $\left\vert {\psi }_{\text{I}}\right\rangle =\left\vert {0}_{1}{1}_{2}{0}_{3}{2}_{4}\right\rangle$ and $\left\vert {\psi }_{\text{O}}\right\rangle =\left\vert {1}_{1}{0}_{2}{1}_{3}{1}_{4}\right\rangle$, then ${S}_{1}^{\text{I}}=0$, ${S}_{2}^{\text{I}}=1$, ${S}_{3}^{\text{I}}=1$, and ${S}_{4}^{\text{I}}=3$, while ${S}_{1}^{\text{O}}=1$, ${S}_{2}^{\text{O}}=1$, ${S}_{3}^{\text{O}}=2$, and ${S}_{4}^{\text{O}}=3$. Then, the matrix ${{\Lambda}}_{\vert {\psi }_{\text{I}}\rangle ,\vert {\psi }_{\text{O}}\rangle }$ reads

$$\begin{equation}{{\Lambda}}_{\vert {\psi }_{\text{I}}\rangle ,\vert {\psi }_{\text{O}}\rangle }=\left(\begin{matrix}\hfill {U}_{1,2}\hfill & \hfill {U}_{3,2}\hfill & \hfill {U}_{4,2}\hfill \\ \hfill {U}_{1,4}\hfill & \hfill {U}_{3,4}\hfill & \hfill {U}_{4,4}\hfill \\ \hfill {U}_{1,4}\hfill & \hfill {U}_{3,4}\hfill & \hfill {U}_{4,4}\hfill \end{matrix}\right).\end{equation} \tag{ 8 }$$

We used a boson-sampling circuit to achieve the full QT of an unknown superposition state in the Fock basis. For obtaining more information of the unknown state, we proposed a simple optical circuit formed by three 50%–50% beam splitters, as shown in figure 1. The unitary matrix U for this circuit is given by

$$\begin{equation}\hat{U}=\left(\begin{matrix}\hfill a\hfill & \hfill -{a}^{2}\hfill & \hfill {a}^{2}\hfill & \hfill 0\hfill \\ \hfill a\hfill & \hfill {a}^{2}\hfill & \hfill -{a}^{2}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {a}^{2}\hfill & \hfill {a}^{2}\hfill & \hfill -a\hfill \\ \hfill 0\hfill & \hfill {a}^{2}\hfill & \hfill {a}^{2}\hfill & \hfill a\hfill \end{matrix}\right)\end{equation} \tag{ 9 }$$

where $a=1/\sqrt{2}$. The input states for the first and the last mode are coherent states given by

$$\begin{equation}\vert {\alpha }_{j}^{(\theta )}\rangle =\sum\limits _{n}^{\infty }\frac{{\text{e}}^{-\frac{1}{2}\vert \alpha {\vert }^{2}}\enspace {\text{e}}^{\text{i}n\theta }\vert \alpha {\vert }^{n}}{\sqrt{n!}}\vert {n}_{j}\rangle ,\end{equation} \tag{ 10 }$$

with θ = 0 for the mode 1 and θ = π/2 for the mode 4. The initial state for mode 2 is the vacuum state $\left\vert {0}_{2}\right\rangle$, and for the mode 3 it is an arbitrary unknown state $\vert {\eta }_{3}\rangle ={\sum }_{\ell =0}^{N}\hspace{2pt}{r}_{\ell }\enspace {\text{e}}^{\text{i}{\phi }_{\ell }}\vert {n}_{3}\rangle$.

Zoom In Zoom Out Reset image size

The considered initial state of our optical circuit is given by

$$\begin{equation}\left\vert {\Psi}\right\rangle =\sum\limits _{m,n=0}^{\infty }\sum\limits _{\ell =0}^{N}{\text{i}}^{n}\frac{{\text{e}}^{-\vert \alpha {\vert }^{2}}\vert \alpha {\vert }^{m+n}}{\sqrt{m!n!}}{r}_{\ell }\enspace {\text{e}}^{\text{i}{\phi }_{\ell }}\left\vert {m}_{1}\right\rangle \left\vert {0}_{2}\right\rangle \left\vert {\ell }_{3}\right\rangle \left\vert {n}_{4}\right\rangle .\end{equation} \tag{ 11 }$$

This state is a superposition of different four-mode states with different number of photons. For the sake of clarity and simplicity, we rewrite the initial state as

$$\begin{equation}\left\vert {\Psi}\right\rangle =\sum\limits _{s=0}^{\infty }\vert {\psi }_{s}^{\text{I}}\rangle ,\end{equation} \tag{ 12 }$$

where $\vert {\psi }_{s}^{\text{I}}\rangle$ is a non-normalized state with the superposition of all the four-mode states with s photons, which reads

$$\begin{equation}\vert {\psi }_{s}^{\text{I}}\rangle =\sum\limits _{\ell =0}^{N}\sum\limits _{n=0}^{s-\ell }{\text{i}}^{n}\frac{{\text{e}}^{-\vert \alpha {\vert }^{2}}\vert \alpha {\vert }^{s-\ell }}{\sqrt{(s-\ell -n)!n!}}{r}_{\ell }\enspace {\text{e}}^{\text{i}{\phi }_{\ell }}\times \left\vert {(s-\ell -n)}_{1}\right\rangle \left\vert {0}_{2}\right\rangle \left\vert {\ell }_{3}\right\rangle \left\vert {n}_{4}\right\rangle .\end{equation} \tag{ 13 }$$

Then, the output state is a superposition of different configurations of how the photons may have arrived to their modes,

$$\begin{equation}\vert {\psi }_{s}^{\text{O}}\rangle =\sum\limits _{{C}_{s}}\hspace{2pt}{\gamma }_{{C}_{s}}\vert ghk{f\rangle }_{{C}_{s}},\end{equation} \tag{ 14 }$$

where C_s denotes a particular configuration with s photons. The probability of measuring a configuration C_s is given by

$$\begin{equation}{P}_{{C}_{s}}={\text{e}}^{-2\vert \alpha {\vert }^{2}}{\left\vert \sum\limits _{\ell =0}^{{N}_{\mathrm{min}}}\sum\limits _{n=0}^{s-\ell }{\text{i}}^{n}\frac{{r}_{\ell }\enspace {\text{e}}^{\text{i}{\phi }_{\ell }}\vert \alpha {\vert }^{s-\ell }\enspace \text{Per}[{{\Lambda}}_{\vert {\psi }_{\text{I}}\rangle ,\vert {\psi }_{\text{O}}\rangle }]}{(s-\ell -n)!n!\sqrt{\ell !g!h!k!f!}}\right\vert }^{2},\end{equation} \tag{ 15 }$$

where N_min = min (s, N). Next, we will consider particular cases of the maximum number of photons of the state $\left\vert {\eta }_{3}\right\rangle$ (N = 1 and N = 2), and then we will extend them to an arbitrary superposition of Fock states.

3.1. The case N = 1

The unknown state is $\vert {\eta }_{3}\rangle ={r}_{0}\vert 0\rangle +{r}_{1}\enspace {\text{e}}^{\text{i}{\phi }_{1}}\vert 1\rangle$ when N = 1. Firstly, for s = 0 we have only one possible output ${C}_{0}=\left\vert 0000\right\rangle$ with probability given by equation (15),

$$\begin{equation}{P}_{{C}_{0}}={\text{e}}^{-2\vert \alpha {\vert }^{2}}{r}_{0}^{2}\end{equation} \tag{ 16 }$$

where we consider ϕ₀ = 0. For s = 1 we have four different outputs, ${C}_{1}\in \left\{\left\vert 1000\right\rangle ,\left\vert 0100\right\rangle ,\left\vert 0010\right\rangle ,\left\vert 0001\right\rangle \right\}$, with the general formula to get the probability given by

$$\begin{equation}{P}_{{C}_{1}}={\text{e}}^{-2\vert \alpha {\vert }^{2}}{\left\vert \vert \alpha \vert {r}_{0}(\text{Per}[{{\Lambda}}_{(\left\vert 1000\right\rangle ,\left\vert {C}_{1}\right\rangle )}]+\text{i}\enspace \text{Per}\left[\right. {\Lambda}\left[\right. (\left\vert 0001\right\rangle \left\vert {C}_{1}\right\rangle ,\left. \right]))+{r}_{1}\enspace {\text{e}}^{\text{i}{\phi }_{1}}\enspace \text{Per}\left[\right. {\Lambda}[(\left\vert 0010\right\rangle ,\left\vert {C}_{1}\right\rangle )]\right\vert }^{2}.\end{equation} \tag{ 17 }$$

The probabilities are

$$\begin{equation}{P}_{\left\vert 1000\right\rangle }=\frac{1}{2}\enspace {\text{e}}^{-2\vert \alpha {\vert }^{2}}\left(\vert \alpha {\vert }^{2}{r}_{0}^{2}+\frac{1}{2}{r}_{1}^{2}+\sqrt{2}\vert \alpha \vert {r}_{0}{r}_{1}\enspace \mathrm{cos}\enspace {\phi }_{1}\right),\end{equation} \tag{ 18 }$$

$$\begin{equation}{P}_{\left\vert 0100\right\rangle }=\frac{1}{2}\enspace {\text{e}}^{-2\vert \alpha {\vert }^{2}}\left(\vert \alpha {\vert }^{2}{r}_{0}^{2}+\frac{1}{2}{r}_{1}^{2}-\sqrt{2}\vert \alpha \vert {r}_{0}{r}_{1}\enspace \mathrm{cos}\enspace {\phi }_{1}\right),\end{equation} \tag{ 19 }$$

$$\begin{equation}{P}_{\left\vert 0010\right\rangle }=\frac{1}{2}\enspace {\text{e}}^{-2\vert \alpha {\vert }^{2}}\left(\vert \alpha {\vert }^{2}{r}_{0}^{2}+\frac{1}{2}{r}_{1}^{2}-\sqrt{2}\vert \alpha \vert {r}_{0}{r}_{1}\enspace \mathrm{sin}\enspace {\phi }_{1}\right),\end{equation} \tag{ 20 }$$

$$\begin{equation}{P}_{\left\vert 0001\right\rangle }=\frac{1}{2}\enspace {\text{e}}^{-2\vert \alpha {\vert }^{2}}\left(\vert \alpha {\vert }^{2}{r}_{0}^{2}+\frac{1}{2}{r}_{1}^{2}+\sqrt{2}\vert \alpha \vert {r}_{0}{r}_{1}\enspace \mathrm{sin}\enspace {\phi }_{1}\right).\end{equation} \tag{ 21 }$$

Note that the parameters r₀, r₁, and ϕ₁ are encoded in the probabilities. Therefore, these parameters can be obtained by calculating the sum or difference of the probability distributions.

3.2. The case N = 2

In this case, we consider $\vert {\eta }_{3}\rangle ={r}_{0}\vert 0\rangle +{r}_{1}\enspace {\text{e}}^{\text{i}{\phi }_{1}}\vert 1\rangle +{r}_{2}\enspace {\text{e}}^{\text{i}{\phi }_{2}}\vert 2\rangle$; the probabilities for s < 2 are the same as in the previous case. For s = 2, we have

$$\begin{align}{P}_{{C}_{2}}& =\frac{{\text{e}}^{-2\vert \alpha {\vert }^{2}}}{g!h!k!f!}\left[\frac{\vert \alpha {\vert }^{2}{r}_{0}}{2}\left(\text{Per}[{{\Lambda}}_{(\left\vert 2000\right\rangle ,\left\vert {C}_{2}\right\rangle )}]-\text{Per}[{{\Lambda}}_{(\left\vert 0002\right\rangle ,\left\vert {C}_{2}\right\rangle )}]+2\text{i}\enspace \text{Per}[{{\Lambda}}_{(\left\vert 1001\right\rangle ,\left\vert {C}_{2}\right\rangle )}]\right)\right.\\ & \quad {\left.+\vert \alpha \vert {r}_{1}\enspace {\text{e}}^{{\imath}{\phi }_{1}}\left(\text{Per}[{{\Lambda}}_{(\left\vert 1010\right\rangle ,\left\vert {C}_{2}\right\rangle )}]+\text{i}\enspace \text{Per}[{{\Lambda}}_{(\left\vert 0011\right\rangle ,\left\vert {C}_{2}\right\rangle )}]+\frac{{r}_{2}\enspace {\text{e}}^{{\phi }_{2}}}{\sqrt{2}}\enspace \text{Per}[{{\Lambda}}_{(\left\vert 0020\right\rangle ,\left\vert {C}_{2}\right\rangle )}]\right)\right]}^{2}.\end{align} \tag{ 22 }$$

For ℓ = 0 and g + h = 2 − n, the permanent is given by

$$\begin{equation}\text{Per}\left[\right. {{\Lambda}}_{\left\vert {(2-n)}_{1}\right\rangle \left\vert {0}_{2}\right\rangle \left\vert {0}_{3}\right\rangle \left\vert {n}_{4}\right\rangle ,\vert {C}_{2}\rangle \left. \right]}={(-1)}^{k}{a}^{2}(2-n)!n!\end{equation} \tag{ 23 }$$

For ℓ ≠ 0, we have the case g + h − d = 2 − ℓ − n and k + f − q = n, where d + q = ℓ. Here, the corresponding permanent reads

$$\begin{align}\text{Per}[{{\Lambda}}_{\left\vert {(2-\ell -n)}_{1}\right\rangle \left\vert {0}_{2}\right\rangle \left\vert {\ell }_{3}\right\rangle \left\vert {n}_{4}\right\rangle ,\vert {C}_{2}\rangle }]& ={B}_{\ell }^{d}\cdot {a}^{2+\ell }\cdot \left[\sum\limits _{x=0}^{d}{B}_{d}^{x}{(-1)}^{x}(2-\ell -n)!\frac{g!}{[g-(d-x)]!}\frac{h!}{(h-x)!}\right]\cdot \\ & \quad \times \left[\sum\limits _{y=0}^{q}{B}_{q}^{y}{(-1)}^{y}n!\frac{f!}{[f-(q-y)]!}\frac{k!}{(k-y)!}\right],\end{align} \tag{ 24 }$$

where ${B}_{j}^{k}=\frac{j!}{k!(j-k)!}$ is the binomial coefficient. For example, for ℓ = 1, $\text{Per}[{{\Lambda}}_{\left\vert 1010\right\rangle ,\left\vert 2000\right\rangle }]$ is given by

$$\begin{equation}\text{Per}[{{\Lambda}}_{\left\vert 1010\right\rangle ,\left\vert 2000\right\rangle }]={(-1)}^{k}{a}^{3}(g-h).\end{equation} \tag{ 25 }$$

When ℓ = 2, ${\text{Per}}_{{{\Lambda}}_{\left\vert 0020\right\rangle ,\left\vert 1001\right\rangle }}$ can be written as

$$\begin{equation}\text{Per}[{{\Lambda}}_{\left\vert 0020\right\rangle ,\left\vert 1001\right\rangle }]=2\cdot {(-1)}^{k}{a}^{4}\left[(f-k)(g-h)\right].\end{equation} \tag{ 26 }$$

Then, the probabilities for the configurations $\left\vert 2000\right\rangle$ and $\left\vert 0020\right\rangle$ are given by

$$\begin{align}{P}_{\left\vert 2000\right\rangle }& =\frac{1}{4}\enspace {\text{e}}^{-2\vert \alpha {\vert }^{2}}\left(\frac{\vert \alpha {\vert }^{4}{r}_{0}^{2}}{2}+\vert \alpha {\vert }^{2}{r}_{1}^{2}+\frac{{r}_{2}^{2}}{4}+\sqrt{2}\vert \alpha {\vert }^{3}{r}_{0}{r}_{1}\enspace \mathrm{cos}({\phi }_{1})\right.\\ & \quad \left.+\frac{\sqrt{2}\vert \alpha {\vert }^{2}{r}_{0}{r}_{2}\enspace \mathrm{cos}({\phi }_{2})}{2}+\vert \alpha \vert {r}_{1}{r}_{2}\enspace \mathrm{cos}({\phi }_{1}-{\phi }_{2})\right),\end{align} \tag{ 27 }$$

$$\begin{align}{P}_{\left\vert 0020\right\rangle }& =\frac{1}{4}\enspace {\text{e}}^{-2\vert \alpha {\vert }^{2}}\left(\frac{\vert \alpha {\vert }^{4}{r}_{0}^{2}}{2}+\vert \alpha {\vert }^{2}{r}_{1}^{2}+\frac{{r}_{2}^{2}}{4}+\sqrt{2}\vert \alpha {\vert }^{3}{r}_{0}{r}_{1}\enspace \mathrm{sin}({\phi }_{1})\right.\\ & \quad \left.-\frac{\sqrt{2}\vert \alpha {\vert }^{2}{r}_{0}{r}_{2}\enspace \mathrm{cos}({\phi }_{2})}{2}-\vert \alpha \vert {r}_{1}{r}_{2}\enspace \mathrm{sin}({\phi }_{1}-{\phi }_{2})\right).\end{align} \tag{ 28 }$$

For the two-photon configuration (s = 2), we can code the different phase differences ϕ_j − ϕ_k and probabilities amplitudes r_j for j, k ∈ {0, 1, 2}.

The same process can be considered for each configuration with a fix number of photons s, and we can find that the unknown phases and amplitudes will be encoded in the different probabilities. This will produce more complex mathematical expressions, thus being hard to calculate the analytical expression to recover some information from the unknown state $\left\vert {\eta }_{3}\right\rangle$. Nevertheless, it looks a suitable problem for ML methods, which can learn and recognize patterns from data. It is important to mention that when the reference states are quantum superpositions $(\left\vert 0\right\rangle +{\text{e}}^{\text{i}\theta }\left\vert 1\right\rangle )/\sqrt{2}$ for θ = 0 and θ = π/2, the mathematical form of the output probabilities are simple. In this case, analytical expressions for all the amplitudes and phases may be found in a straightforward manner, reaching full tomography with 100% accuracy.

Entanglement is one of the most important resources for quantum information. It describes nonlocal correlations between quantum states, and it has become an important tool for understanding the states of many-body systems. For bipartite systems, the entanglement entropy has become a theoretical measure for categorizing such states. For the entanglement coding, we used the optical circuit shown in figure 2, where there is an unknown two-mode input state. We will make use of the von Neumann entropy of a state, which is defined as

$$\begin{equation}S(\rho )=-\mathrm{Tr}(\rho \enspace \mathrm{ln}\enspace \rho ),\end{equation} \tag{ 29 }$$

where ρ is the density operator of the composite system, while Tr(⋅) denotes the trace.

Zoom In Zoom Out Reset image size

For the calculation of the initial state, we suppose that the first and last modes are coherent states given by

$$\begin{equation}\vert {\alpha }_{j}^{(\theta )}\rangle =\sum\limits _{n=0}^{\infty }\enspace \frac{{\text{e}}^{-\frac{1}{2}\vert \alpha {\vert }^{2}}{\text{e}}^{\text{i}n\theta }\vert \alpha {\vert }^{n}}{\sqrt{n!}}\vert {n}_{j}\rangle ,\end{equation} \tag{ 30 }$$

with θ = 0 for the mode 1, θ = π/2 for the mode 4 as in the previous case, and a bipartite system for the second and third mode, it means a wave function of the following form

$$\begin{equation}\vert \psi \rangle =\sum\limits _{j,v=0}^{N}{r}_{jv}\enspace {\text{e}}^{\text{i}{\phi }_{jv}}\vert {j}_{2}\rangle \vert {v}_{3}\rangle ,\end{equation} \tag{ 31 }$$

where the parameters r_jv and ϕ_jv satisfy ${\sum }_{j,v=0}^{N}\hspace{2pt}{r}_{jv}^{2}=1$ and ϕ_jv ∈ [0, 2π]. Then, the initial state of our optical circuit is

$$\begin{equation}\left\vert {\Psi}\right\rangle =\sum\limits _{m,n=0}^{\infty }\sum\limits _{j,v=0}^{N}{\text{i}}^{n}\frac{{\text{e}}^{-\vert \alpha {\vert }^{2}}\vert \alpha {\vert }^{m+n}}{\sqrt{m!n!}}{r}_{jv}\enspace {\text{e}}^{\text{i}{\phi }_{jv}}\left\vert {m}_{1}\right\rangle \left\vert {j}_{2}\right\rangle \left\vert {v}_{3}\right\rangle \left\vert {n}_{4}\right\rangle .\end{equation} \tag{ 32 }$$

We describe the bipartite system in modes 3 and 4 by the density operator ρ_AB = |ψ⟩⟨ψ|, and the reduced density matrix of the subsystem A(B) as ρ_A = Tr_B (ρ_AB ) (ρ_B = Tr_A (ρ_AB )). ${\mathrm{Tr}}_{\mathcal{S}}$ means tracing over subsystem $\mathcal{S}$. The entanglement E(ρ_AB ) of a bipartite system ρ_AB may be defined as the von Neumann entropy of the reduced density matrix of a subsystem,

$$\begin{equation}E({\rho }_{AB})=S({\rho }_{A})=S({\rho }_{B}).\end{equation} \tag{ 33 }$$

For simplicity, we rewrite the state given by equation (32) as

$$\begin{equation}\left\vert {\Psi}\right\rangle =\sum\limits _{s=0}^{\infty }\vert {\psi }_{s}^{\text{I}}\rangle ,\end{equation} \tag{ 34 }$$

where $\vert {\psi }_{s}^{\text{I}}\rangle$ is a non-normalized state of the superposition of all the four-mode states with s photons, which reads

$$\begin{equation}\vert {\psi }_{s}^{\text{I}}\rangle =\sum\limits _{j,v=0}^{N}\sum\limits _{n=0}^{s-j-v}{\text{i}}^{n}\frac{{e}^{-\vert \alpha {\vert }^{2}}\vert \alpha {\vert }^{s-j-v}}{\sqrt{(s-j-v-n)!n!}}{r}_{jv}\enspace {\text{e}}^{\text{i}{\phi }_{jv}}\left\vert {(s-j-v-n)}_{1}\right\rangle \left\vert {j}_{2}\right\rangle \left\vert {v}_{3}\right\rangle \left\vert {n}_{4}\right\rangle .\end{equation} \tag{ 35 }$$

Then, the output state is a superposition of the different configurations of how photons may have arrived to output modes,

$$\begin{equation}\vert {\psi }_{s}^{\text{O}}\rangle =\sum\limits _{{C}_{s}}\hspace{2pt}{\gamma }_{{C}_{s}}\vert ghk{f\rangle }_{{C}_{s}},\end{equation} \tag{ 36 }$$

where C_s denotes a particular configuration with s photons. The probability of measuring a specific configuration C_s reads

$$\begin{equation}{P}_{{C}_{s}}={\text{e}}^{2\vert \alpha {\vert }^{2}}{\left\vert \sum\limits _{j,v=0}^{N}\sum\limits _{n=0}^{s-j-v}\frac{{i}^{n}\vert \alpha {\vert }^{s-j-v}\enspace {\text{e}}^{\text{i}{\phi }_{jv}}{r}_{jv}\enspace \text{Per}[{{\Lambda}}_{\left\vert {s}^{\ast }jvn\right\rangle ,\left\vert ghkf\right\rangle }]}{(s-j-v-n)!n!\sqrt{g!h!k!f!j!v!}}\right\vert }^{2},\end{equation} \tag{ 37 }$$

where s* = s − j − v − n. Now, each set of probabilities is related to a given degree of entanglement. Then, again via the use of an ML protocol based on training a pattern recognition algorithm, we can estimate this essential feature.

ML was used to find the relation between the patterns and the initial arbitrary state so that the model, once trained, could infer the state of a new pattern not belonging to the training set. In particular, we developed regression models for the values of amplitude and phase, calculating the associated fidelity as figure-of-merit. Due to the sparsity of the data set, support vector regressors (SVRs) were chosen to carry out the regression of amplitudes and phases as our first approach [43, 44]. SVRs are a generalization of support vector machines, which try to solve classification problems by formulating them as convex optimization problems. In this manner, one has to find the suitable hyperplanes that classify correctly as many training samples as possible. In particular, SVRs work by creating a transformed data space in which the problem is more easily solvable and, ideally the problem is transformed into a linear one. That transformation between spaces is carried out by the so-called kernels. Gaussian, linear, and polynomial kernels have been used in this experimentation. An SVR introduces a region in the hyperspace of the problem called -tube, within which, all predictions are considered as correct.

The first three tables show the results of the predictions made by the ML models, evaluating those predictions by means of the fidelity as figure-of-merit. For each data instance, there is a corresponding fidelity value; the tables show the mean value of the fidelities as well as its standard deviation and the quartiles as measures of dispersal. A percentile X% means that in X% of our data samples, fidelities are below the percentile X% value. For example, in table 1, for the Radial basis function (RBF) (Gaussian SVR), 25% of the samples give fidelities below 0.6397; this of course means that 75% of the results give fidelity values above 0.6397. Tables 4 and 5 report the entanglement estimation. The entropy mean value is given for reference purposes to compare its magnitude and its mean absolute error (MAE). We also give the standard deviation of the absolute errors as well as their quartiles and the coefficient of determination (R²-score) for the models. The coefficient of determination evaluates the quality of the regression, being 1 its best possible value; a constant model that gives the same output value disregarding the input features would have a score of 0.

Starting with the specific results of each experiment, table 1 shows the values of fidelity obtained by SVRs implemented with Gaussian (RBF), linear, and polynomial kernels for the case of N = 3. The maximum fidelities are achieved by the Gaussian kernel, that will be hence selected for ulterior analyses. The polynomial kernel obtains slightly lower fidelities than the Gaussian one, whereas the linear kernel shows the poorest modeling capability among the three.

For the sake of comparison with SVRs, we considered a state-of-the-art method, such as extremely randomized trees (ERTs) [45], which build multiple regression trees. Each tree takes a random subset of the input features, while nodes are randomly split for the whole data set (in contrast with the well-known random forest, there is no bootstrap). To reduce the problem dimensionality, we made use of a principal component analysis (PCA) [37]. We set the explained variance ratio at 0.999, so that we could get a good-enough representation of our data that allows for a reliable fidelity modeling while considerably reducing the computational time.

Table 2 shows the fidelity values obtained by SVRs with a Gaussian kernel, for non-PCA and PCA versions of the data sets. The use of PCA seems to have no visible effects on the fidelity performance. In fact, fidelities tend to be slightly higher when using PCA, thus suggesting its suitability for dimensionality reduction. Table 2 also shows that the most reliable fidelities are achieved for N = 2, slightly lower values are obtained for N = 1, and even lower for N = 3.

Table 3 includes the values of fidelity obtained by ERTs for both non-PCA and PCA versions of the data sets. As in the case of SVRs, the use of PCA does not affect the performances in a great extent. In particular, the use of PCA provides slightly better results for N = 2 and N = 3, while for N = 1 the use of the original data without PCA yields higher values of fidelity. The comparison of tables 2 and 3 show that the ERTs carry out a better regression than SVRs for N = 2 and N = 3, being an SVR approach better for N = 1. This might likely be due to the fact that N = 1 corresponds to a more sparse data set, where the SVR is a more adequate and natural solution. Furthermore, that data set may not contain enough variability to exploit ERT capabilities. In summary, the use of SVRs is suggested for N = 1 and ERTs for any N > 1. Furthermore, as PCA does not have a relevant impact on performance, tends to lead to higher fidelities, and allows a reduction of the computational times, its use to reduce the dimensionality of the data sets is encouraged.

Results about entanglement estimation via the corresponding entropy are shown in tables 4 and 5. The performance of SVRs and ERTs are compared via the MAE and R² scores for different data sets. ERTs outperform SVRs in all cases except for N = 1, when there is no PCA preprocessing again very likely due to sparsity. This conjecture is reinforced by the fact that after using PCA, hence reducing sparsity, SVRs are not better regressors than ERTs, even for N = 1. As in the estimation of fidelities, PCA does not have a relevant impact on performance. In fact, its use tends to lead to slightly better performances (higher R²-scores and lower MAEs) whilst allowing a reduction of the computational times. Therefore, its use to reduce the dimensionality of the data sets is encouraged. We may infer that, for entanglement estimation, ERTs represent a more adequate choice than SVRs.

We have proposed a state characterization via a pattern recognition algorithm in ML for photonics. It is based on the estimation of quantum features by using the output photon number distribution in a photonic circuit, similar to boson-sampling protocols for continuous variables. We have focused on the estimation of single-mode phases and amplitudes and, also, on the two-mode entanglement estimation.

SVRs and ERTs have been used to extrapolate the estimation for states not present in the training set, i.e. new probability patterns. The obtained fidelities hint that ML estimations are reliable and can be used to boost the proposed QT protocol. In particular, the use of ERTs with previous PCA preprocessing seems to be a suitable approach for setups with N > 1.

The same two ML approaches were employed for entanglement estimation through modeling the von Neumann entropy, yielding R² scores higher than 0.75, thus suggesting the appropriateness of ML for this relevant task. In this case, the use of ERTs with a previous PCA preprocessing turns out to be the most adequate choice for all cases.

The proposed photonic circuit with four modes and three beam splitters is experimentally accessible. It means that by using a bigger circuit, with more output probabilities for the number of photons, ML algorithms will likely increase its performance, since the patterns will be more complex and will encode more information. Nevertheless, the numerical generation of the data will be a hard task due to the complexity of the problem. Also, the generation of experimental data for the training set is a demanding task due to the difficulties to characterize a quantum state. In any case, it is a useful technique even if we use simple photonic circuits or historical data produced by current photonic experiments for more complex setups. Finally, this work shows that ML techniques can be suitably used for state characterization without the burden of full tomography, paving the way for more sophisticated tools that may help for fast estimation of quantum features.

We acknowledge support from QMiCS (820505) and OpenSuperQ (820363) Projects from EU Flagship on Quantum Technologies, National Natural Science Foundation of China Grant (NSFC) (12075145), Shanghai Government Grant STCSM (2019SHZDZX01-ZX04), Spanish Government Grant PGC2018-095113-B-I00 (MCIU/AEI/FEDER, UE), Basque Government IT986-16, EU FET Open Quromorphic and EPIQUS Projects.

The data that support the findings of this study are available upon reasonable request from the authors.