1. Introduction

The use of classical and quantum properties of light opens up new possibilities for information processing methods [1]. Programmable optical circuits play an important role in classical and quantum optical communications [2–5], in quantum computing [6–8], and in machine learning problems [9–11]. The key element of such devices is a programmable linear optical interferometer. When control parameters are fixed, it performs a linear transformation of $N$ input fields:

One can tune the interferometer to some target unitary transformation $V$ in two ways. The first way is empirical. For each $V$, one optimizes the control parameters in order to bring the experimental results of optical fields measurements closer to the target ones [17]. This does not require building an exact model of interferometer internal structure, but is impractical, since it is necessary to measure all matrix elements at each step. For each new target matrix $V$, the procedure must be repeated. Another way is to first build a complete mathematical model of the interferometer based on the results of some set of measurements [18]. Then, the search for the desired set of control parameters is performed without any additional measurements. This procedure faces some difficulties. Since large-scale linear optical devices can be implemented in the integral design only, it is almost impossible to split the large-scale reconstruction problem into several low-scale ones. It needs additional optical outputs or modular optical circuits [19] that causes additional losses and instability. Trying to characterize an integrated device as a whole, one has to face the fact that the results of each measurement are determined by a large number of internal device parameters. As a result, the efficiency of the optimization procedure rapidly decreases with increasing dimension $N$.

In this paper, we present a method for constructing a model of a programmable integrated interferometer based on algebraic operations only. This approach needs limited number of operations and is much more computationally efficient than the one using numerical optimization. The method is based on the layered interferometer architecture [15], which is described in Section 2. Further in Section 3, we describe the procedure for reconstructing interferometer mixing layers. As in [18], it requires the measurement of all matrix elements of the complex transfer matrix $V$ (up to a global phase) for a given configuration of phase layers. Common methods for measuring the transfer matrices of optical interferometers using laser [20], thermal [21], and two-photon [22,23] fields make it possible to measure matrix elements up to input and output phases only. A non-ambiguous measurement of the matrix $V$ can be performed using homodyne detection of all the input and output coherent fields in Eq. (1) and requires $O(N^2)$ measurements [24]. In Section 4, we show that the proposed algebraic approach allows one to predict the transformation matrix $V$ for an arbitrary configuration of phase layers with an accuracy not lower than that for the method based on numerical optimization. The mixing layers are determined ambiguously, however, in the case of homogeneous losses, our method makes it possible to uniquely determine the normalized transmission coefficients of each individual layer. In addition, we show that the method is highly robust to phase layers control errors.

2. Interferometer model

Below we consider the universal layered interferometer architecture consisting of alternating mixing and phase layers [15,18] (Fig. 1). The corresponding transfer matrix has the form

Here $U_k$ and $\Phi _k$ are transfer matrices of the mixing and phase layers, respectively. Each mixing layer is an $N$-channel transformation that sufficiently distributes the signal of each input channel to all the output channels. All matrices $U_k$ can be different. Ideally, they are unitary, but this property is lost in the presence of linear losses. For the phase layer, each channel is subjected to a controlled phase shift through a thermo- or electro-optical effect. The corresponding transfer matrix has a diagonal form:

Expression (2) can also be used to describe any other interferometer architecture (for example, the Reck [12] or Clements [13] architecture). In this case, $U_k$ can be the transfer matrices of beam splitters acting on two channels only, and the phase layers do not control all channels, but only a few ones.

 

Fig. 1. Schematic representation of the layered architecture of a programmable interferometer. The number of mixing layers $K$ is equal to the number of channels $N=5$.

Download Full Size | PDF

With the number of mixing layers $K=N$, the layered architecture makes it possible to obtain an arbitrary transformation $V$ due to the proper choice of phases $\varphi _n^{(k)}$. We mostly consider this particular case.

The problem of constructing an interferometer model is reduced to determining all transfer matrices $U_k$ and $\Phi _k$. The dependence of $\varphi _n^{(k)}$ on the control parameters can be pre-calibrated even for unknown $U_k$. For this, a simple $\sin \varphi$ dependence of the interferometer output intensity on the value $\varphi$ of any single controlled phase is used (all other phases are disabled) [2]. Below we assume that all phase layers are pre-calibrated (later we will show that our reconstruction method is robust to phase control errors). Reconstruction of the matrices $U_k$ is complicated by the fact that there is no way to “turn off” all mixing layers, except for one, and perform its measurement. However, in the next section, we show that it is possible to obtain a matrix equation that depends on a single mixing layer only by an appropriate choice of phase layers configurations. This equation is reduced to an eigendecomposition problem, which is solved algebraically.

3. Mixing layers reconstruction

In this section, we assume that we have full control over the phase layers. Later in Section 4, we show that the proposed method is robust to control errors. We also assume that none of the mixing layers has absolute losses in any channel, so all the matrices $U_1, \dots U_K$ are invertible. We also do not restrict them to be unitary.

We describe the reconstructing algorithm in several steps. First, we consider the ideal case when, for any configuration of the interferometer, we know exactly its full transfer matrix $V$. Next, we consequentially take into account the impossibility of experimental determination of the global phase of $V$, the presence of errors in measurements and phase layers control, and the special case of homogeneous losses inside the interferometer. Finally, we present the pseudo-code for the complete algorithm. Its software implementation is available online [25].

3.1 Ideal case

Consider two configurations of the interferometer phase layers. The first configuration has all the phase layers disabled. The corresponding transfer matrix is $V_0=U_K\dots U_2U_1$. Next, the phases $\{\varphi _n, n=1,\dots,N\}$ are set in the second phase layer: $V_1=U_K\dots U_2\Phi U_1$, where $\Phi =\operatorname {diag}[\exp (i\varphi _n)]$. Consider the matrix $A_1=V_0^{-1}V_1=U_1^{-1}\Phi U_1$. Since the matrix $\Phi$ is known and diagonal, this equation is solved with respect to $U_1$ by the diagonalization of $A_1$. By induction, one can iteratively obtain all mixing matrices up to $U_{K-1}$. In particular, we enable only the $(k+1)$-th phase layer to obtain the $k$-th mixing layer:

Using $V_k$ and the results for all previous $k-1$ matrices we compute the following matrix:

The following equation determines its eigendecomposition:

Since eigenvectors diagonalize $A_k$, we can use it to construct $U_k^{-1}$. To do this, one must sort the solutions of Eq. (6) in such a way that matches eigenvalues with the diagonal of $\Phi$. Let us define the sorting indices $S=\{s_n: \lambda _{s_n}=\exp (i\varphi _n)\}$. Then the inverse transfer matrix has the form

This procedure goes successively for $k=1,\dots,K-1$. Finally, the transfer matrix of the last mixing layer is $U_K=V_0U_1^{-1}U_2^{-1}\dots U_{K-1}^{-1}$. The procedure does not involve the first and last phase layers.

Note that for the correct work of the described method it is necessary that all $\lambda _q$ from Eq. (6) (so the phases $\varphi _n$) must be different. Eigenvectors with equal eigenvalues cannot be uniquely defined, this leads to errors in constructing the columns of $U_k^{-1}$.

3.2 Ambiguity

The described algorithm determines all mixing layers of the interferometer up to following transformations:

 

Fig. 2. Illustration of the impossibility to distinguish the losses and phase shifts at the output of one mixing layer from the losses and phase shifts at the input of the next mixing layer. Boxes with symbols $\theta$ and $\eta$ denote phase shifts and linear losses, respectively.

Download Full Size | PDF

3.3 Effect of the global phase

The approach described above gives an exact solution, when the matrices $V_k$ are exactly known. In practice, the elements of the transfer matrices are measured experimentally and, therefore, differ from the original ones. These differences are primarily related to the undetectable global phase. In this regard, in Eq. (6), we obtain eigenvalues $\lambda _q=\exp (i\delta _q)$, where the set of phases $\{\delta _q\}$ differs from the original set $\{\varphi _n\}$ by some constant value $\chi$. It makes the matching between $\{\delta _q\}$ and $\{\varphi _n\}$ ambiguous and could violate the sorting of eigenvectors. For example, a set of phases

One can observe that this set is invariant under the full $2\pi$ rotations only. Further, it is convenient to move from absolute phases to relative phases between adjacent channels of the phase layer (channel $N$ is considered adjacent to channel $1$):

It is easy to see that the relative phase for the first channel is twice the others for the choice (Eq. (11)):

 

Fig. 3. The distribution of phases on the unit circle. (a) Uniform distribution of controllable phases on the single phase layer for $N = 5$. If the constant $\chi$ is a multiple of $2\pi /5$ then the transformation $\varphi _n \rightarrow \varphi _n+\chi \ (\mathrm {mod}\ 2\pi )$ leads to a similar distribution making it impossible to unambiguously compare the eigenvalues $\lambda _q=\exp (i\delta _q)$ with the phases $\varphi _n$. (b) The choice of asymmetric phase distribution (11) distinguishes $\varphi _1$ from the others by the greatest distance to the previous phase. (c) The search for the largest distance between adjacent phases of eigenvalues $\lambda _q$ allows one to determine the index $q$ corresponding to the phase $\varphi _1$ and retrieve the constant $\chi$. (d) Strong phase fluctuations of the eigenvalues $\lambda _q$ at large $N$ (in this case $N = 19$) violates the sorting of the eigenvectors.

Download Full Size | PDF

Systematic and statistical errors also distort the spectrum of the matrix $A_k$, leading to more general eigenvalues: $\lambda _q=r_q\exp (i\delta _q)$, where $r_q\approx 1$, $\delta _q\approx \varphi _q+\chi$. If the errors are not too large then we can again use the above approach, but instead of an exact equality we look for an approximate one:

Note that in practice this minimization is not necessary. The set Eq. (11) is sorted in ascending order, so one could also sort all $\delta _q$ and then find the maximum value of $\Delta \delta _q$. In view of Eq. (13), this value corresponds to the index $s_1$. The remaining elements will follow it on the unit circle in the anti-clockwise direction (Fig. 3(c)).

3.4 Noise handling

The described algorithm works well on small dimensions. However, as $N$ increases, the density of $\{\varphi _n\}$ values on the unit circle increases. As a result, the eigenvalue sorting rule becomes more noise sensitive (Fig. 3(d)). One can reduce this effect by activating only a few phases at a time and calculating only a few columns of the matrix $U_k^{-1}$. Denote by $C=\{c_m\}$ the set of indices of these columns and select the following set of phases (Fig. 4(a)):

 

Fig. 4. An example of a partial reconstruction of the inverse matrix of the $k$-th mixing layer. The dimension of the interferometer is $N=9$. $M = 4$ matrix columns from 5 to 8 are estimated. (a) The $k$-th phase layer configuration. (b) The phases of matrix $A_k$ eigenvalues are sorted on the unit circle in ascending order. The first element is considered to be the one that is separated from the previous one by the largest distance. The first $N-M=5$ elements correspond to eigenvectors with indices $n\notin C$. The remaining $M = 4$ elements are eigenvectors with indices $n\in C$. (c) The corresponding eigenvectors are set to be the columns of the matrix $U_k^{-1}$.

Download Full Size | PDF

The smaller the number $M$ of simultaneously estimated columns of $U_k^{-1}$, the less densely the corresponding phase values can be located. As a result it is possible to resolve the individual eigenvalues of the matrix $A_k$ from Eq. (6) quite well even for large values of dimension $N$. In the limit $M=1$, we reconstruct the columns of the matrix $U_k^{-1}$ one by one. This, however, complicates the experiment: it is required to measure $\propto N/M$ transfer matrices for each mixing layer. The total number of transfer matrices used is $L=(K-1)\lceil N/M \rceil + 1$.

Note that in the above algorithm, the phases $\delta _q$ of $A_k$ eigenvalues were used for the correct eigenvectors sorting only. At the same time, the absolute values of $\delta _q$ contain information about the true values of the controlled phases $\varphi _n$ at the phase layer. These values can serve for more accurate calibration of the phase layer, taking its cross-talk into account.

3.5 Case of homogeneous losses

The described procedure does not impose any constraints on resulting matrices $U_k$. In certain cases, it may be known a priory that the interferometer losses are homogeneous. Then each $U_k$ belongs to the set of unitary matrices up to a constant factor $\eta _0$ describing the homogeneous losses. Thus, one can enhance the accuracy by projecting each $U_k$ to this set. For some non-unitary matrix $X$, the projection is given by

Here matrices $Q$, $S$ and $W$ form the singular value decomposition of $X$. The matrices $Q$ and $W$ are unitary, while $S$ is non-negative diagonal matrix and equals to the identity matrix for unitary $X$. Within the unitary model, the transformations $\Lambda _k$ from Eq. (8) are limited by arbitrary phase shifts only (there are no linear loss blocks in Fig. 2). This ambiguity does not affect the normalized transmission coefficients of mixing layers

Thus, this coefficients are defined unambiguously. Below we consider the case $\eta _0=1$ since homogeneous losses could be measured separately and compensated by collecting more data.

3.6 Algorithm pseudo-code

Algorithm 1 summarizes the described procedure of mixing layers reconstructions. It takes the following parameters as the input: $N$ – the number of input and output channels of the interferometer; $K$ – the number of mixing layers; $M$ – the maximal number of $U_k^{-1}$ columns, estimated at once. The output of the algorithm are transfer matrices $U_k$ ($k=1,\dots,K$) of all mixing layers.

The algorithm software implementation is available online [25].

4. Numerical simulation

In this section, we analyze the efficiency of the proposed algorithm in the presence of various kinds of errors. In addition, we compare it with the optimization algorithm from [18].

4.1 Error models

We define random mixing layers in two steps. In the first step, we generate a random unitary matrix

Further we assume $\gamma = 0.1$. In the second step, if necessary, we introduce heterogeneous linear losses:

In other words, we insert random linear losses of no more than $\beta$ into the “middle” part of each random unitary transformation. Random matrices of the form (21) are generated independently for each mixing layer to create a true model of the programmable interferometer. We denote these matrices as $U_k^{\text {true}}$.

Algorithm 1. Mixing layers reconstruction

View Table

We also simulate the errors of phase layers control. Let $\{\varphi _n\}$ be the desired set of phases in the $k$-th phase layer. The relative control error of each phase is given by $\alpha$:

We simulate the statistical noise of the transfer matrix measurements by a random Gaussian error for each matrix element:

4.2 Benchmarks

After the numerical experimental data generation we reconstruct the interferometer model. Then we estimate the ability of the model to predict the full transfer matrix given arbitrary phase layer configuration (each phase of each phase layer is chosen randomly from $0$ to $2\pi$) by calculating the fidelity between the predicted $V_{\text {pred}}$ and true $V_{\text {true}}$ transfer matrices of the entire interferometer (the normalized squared dot product of the matrices):

The quantity $\Delta F = 1 - F$, in turn, plays the role of the distance between the matrices, which is insensitive to the value of the undetectable global phase (similar to the distance between two quantum processes [26]). We consider the $95$-th percentile $\Delta F_{95}$ in a test set of $1000$ random phase layer configurations [27]. Note that the value of $\Delta F_{95}$ is affected not only by the quality of the mixing layers reconstruction, but also by the presence of phase layer control errors in $V_{\text {true}}$, which are assumed to be zero in $V_{\text {pred}}$.

In Section 3.5, we pointed out that the method makes it possible to uniquely determine the normalized transmission coefficients (Eq. (18)) of mixing layers. Therefore, we also consider the measure

It characterizes the maximal error of transmission coefficients estimation over the all mixing layers. Here $T_k^{\text {true}}$ and $T_k^{\text {rec}}$ are respectively the true and reconstructed transmission matrices of the $k$-th mixing layer. Both $\Delta F_{95}$ and $\Delta T_{\max }$ are insensitive to homogeneous losses. Considering $\Delta T_{\max }$ in addition to $\Delta F_{95}$ is useful for two reasons. First, $\Delta F_{95}$ is an integral estimate that includes both errors in determining mixing layers and errors in controlling phase layers. Since the latter remains unknown, $\Delta F_{95}$ may be quite high even if each mixing layer is determined accurately. Secondly, the determination of the transmission coefficient is of independent interest for debugging the manufacturing technology of integrated interferometers. The use of $\Delta T_{\max }$, however, makes sense for the unitary model only (see Section 3.5). In the presence of heterogeneous losses, the matrices $T_k^{\text {rec}}$ are determined ambiguously and it makes no sense to compare them directly with $T_k^{\text {true}}$.

4.3 Comparison with optimization method

Let us consider the case when the number of interferometer channels $N$ and the number of mixing layers $K$ are 5 ($N = K = 5$), and there are no linear losses ($\beta =0$) and phase control errors ($\alpha =w=0$). We generate a random programmable interferometer according to Section 4.1 and perform the reconstruction of all its mixing layers according to Section 3. As in [18], we use the model of unitary mixing layers. The initial data used in the algorithm is also fed to the input of the algorithm from [18] with an open source code. We use the version of this algorithm that uses the approximate prior knowledge of mixing layers transfer matrices. Following the example from the library files, the following algorithm parameters are chosen: mini-batch size is five, the prior knowledge coefficient $\beta =0.24$. The number of training epochs varies from $20$ to $1500$ depending on the dimension $N$ and the measurement protocol (the value is high enough to reliably rich the accuracy saturation).

In total, we considered 3 different measurement protocols:

Figure 5 shows the simulation results. In the absence of phase layers control errors, the dependencies of $\Delta T_{\max }$ (Fig. 5(a)) and $\Delta F_{95}$ (Fig. 5(b)) on the measurement error $\varepsilon$ turned out to be close, regardless of the measurement protocol and the reconstruction method used. The closeness is so high that plotting them in the same axes on a log scale make the distinction barely seen. Approximately, we have observed $\Delta T_{\max }\propto \varepsilon$ and $\Delta F_{95}\propto \varepsilon ^2$. From Fig. 5(c) it can be seen that the presence of phase layer control errors $\alpha$ does not affect our method’s ability to estimate the mixing layer transmission coefficients $[T_k]_{mn}$: it is limited by the interferometer measurements accuracy only. This is because the algorithm uses the phases values only to do the correct sorting. However, starting from a certain critical value of the error $\alpha _c$, the accuracy of determining $[T_k]_{mn}$ sharply decreases, since the sorting is violated. The value of $\alpha _c$ can be increased by decreasing the number of unique phases at each iteration (it is defined by the parameter $M$ of the algorithm). The optimization algorithm from [18] relies on the particular phase values, so the control errors impair the estimation of $[T_k]_{mn}$. Despite this, the algorithm manages to train the entire interferometer as a whole with an error $\Delta F_{95}$ not higher than that provided by our algorithm (Fig. 5(d)). Approximately, we have observed $\Delta F_{95}\propto \alpha ^2$.

 

Fig. 5. The errors in determining the transmission coefficients $\Delta T_{\max }$ (a,c) and predicting the full transfer matrices $\Delta F_{95}$ (b,d). The average values over 100 independent trials are shown. Dimension $N=5$, number of mixing layers $K=5$, no losses ($\beta =0$). The results for our algorithm (solid lines) and for the algorithm from [18] (diamonds) are shown. (a,b) Dependence on the measurement error $\varepsilon$; no phase control errors ($\alpha =w=0$). (c,d) Dependence on the phase control error $\alpha$; measurement error $\varepsilon =10^{-3}$, cross-talk strength $w=0.1$.

Download Full Size | PDF

Despite the closeness of the results for optimization algorithm and our algorithm, the latter performs much faster. To show this, consider two extreme versions of our algorithm ($M = N$ and $M = 1$). For the optimization algorithm from [18], we use a set of 300 training phase layer configurations. Starting from a certain value, changing this number does not affect the computation speed significantly. The number of epochs is chosen to be minimal in order to guarantee a correct result. On the range of considered dimensions $N$, the presented reconstruction algorithm worked more than two orders of magnitude faster (Fig. 6). We also considered the optimization algorithm in the “black box” model (the absence of any prior knowledge about the mixing layers). In this case, as expected, there was a rapid growth in computational complexity with an increase in $N$ from two to six (as in [18], we also failed to obtain the convergence of the optimization algorithm in the “black box” model for $N > 6$). Our algorithm, also based on the “black box” model, showed a much slower increase in computation time. The independent trials of numerical experiments were performed in parallel on 52 CPU cluster; Intel Xeon CPU E5450 @ 3.00GHz.

 

Fig. 6. Comparison of the computation time for different methods versus the dimension $N$. The average values over 100 independent trials are shown. The number of mixing layers $K = N$, no losses ($\beta = 0$), no phase control errors ($\alpha =w=0$), transfer matrices measurement error $\varepsilon =10^{-5}$.

Download Full Size | PDF

4.4 Algorithm scalability

Above we have mentioned that the total number of required interferometer configurations for our method is $L=(K-1)\lceil N/M\rceil +1$, where $N$ is the number of input and output channels, $K$ is the number of mixing layers, and $M\in [1;N]$ is the number of simultaneously estimated columns of the inverse matrix of a single mixing layer. The time required to measure the transfer matrix for each configuration is estimated as $O(N^2)$ [24]. Therefore, for $K=N$ the total measurement time scales as $O(N^4/M)$.

The reconstruction algorithm requires $O(L)$ algebraic operations including matrix multiplication, taking the inverse matrix, and computing matrix eigendecomposition. The complexity of these procedures does not exceed $O(N^3)$ [28]. For $K=N$ the total complexity of the entire procedure can be estimated as $O(N^5/M)$.

From these estimates, one can see the number $M$ to be an important parameter that affects the algorithm speed. The value $M = N$ is optimal from the point of view of measurement and computation speed. However, in this case, with increasing $N$, the sensitivity to control errors of phase layers increases: if the error $\alpha$ is higher than the critical value $\alpha _c$, the mixing layers estimates accuracy sharply decreases. At the same time, for a fixed $M$, the value of $\alpha _c$ changes slightly with increasing $N$ (Fig. 7(a)). Note also that for $\alpha <\alpha _c$ transmission coefficients matrices $T_k$ are determined with approximately the same accuracy for all mixing layers (Fig. 7(b)). That is, despite the iterative nature of the algorithm (see Eq. (5)), there is no effect of error accumulation.

 

Fig. 7. (a) Scalability of the errors $\Delta T_{\max }$ in determining the mixing layers transmission coefficients. The average values over 100 independent trials are shown. Number of mixing layers $K = N$, no losses ($\beta = 0$), cross-talk strength $w=0.1$, transfer matrices measurement error $\varepsilon =10^{-5}$. (b) Distribution of transmission coefficients errors for various mixing layers within a single trial; confidence intervals show the upper and lower quartiles over all $m,n=1,\dots,N$. Parameters: $N=K=25$, $\beta =0$, $w=0.1$, $\alpha =10^{-3}$, $\varepsilon =10^{-5}$. (c) Scalability of the errors $\Delta F_{95}$ in predicting the full transfer matrices. The average values over 100 independent trials are shown. Parameters: $K = N$, $M = N$, $\alpha =w=0$.

Download Full Size | PDF

The prediction error $\Delta F_{95}$ grows polynomially with increasing $N$ (Fig. 7(c)). This is primarily because at a fixed data acquisition time, the size $L$ of the measurement protocol increases. As a result, the measurement accuracy of each transfer matrix decreases (see Eq. (24)), impairing the mixing layers estimation accuracy. The figure shows that one can slightly improve the accuracy by performing the projection of the mixing layers to the set of unitary matrices (Section 3.5). This, however, does not affect the rate of growth (about $\Delta F_{95} \propto N^4$). When the heterogeneous linear losses are present, the rate of growth becomes higher (about $\Delta F_{95} \propto N^{5.8}$ for $\beta =0.1$). This can be explained by the fact that the condition number of the transfer matrices of the mixing layers increases when heterogeneous losses occur. As a result, the matrix inversion procedure becomes more sensitive to statistical fluctuations.

5. Conclusion

We have shown that the estimation of mixing layers in the layered architecture of the programmable integrated optical interferometer is possible without any use of numerical optimization. The proposed algorithm has a single parameter $M$ "- the number of simultaneously estimated columns of a mixing layer inverse transfer matrix. We have performed the numerical simulation and compared the results with the results obtained using the optimization algorithm. It has been found that under the equal conditions both methods provide the same accuracy, regardless of the measurement protocol. At the same time, the exclusively algebraic approach in our algorithm has shown a significantly lower computation time.

Moreover, our algorithm uses the control phases values only to correctly sort the columns of the adjacent mixing layer. This makes the method almost insensitive to phase control errors up to some critical error value. This critical value can be increased by choosing a smaller value of the parameter $M$. Small $M$, however, increases the number of operations in the reconstruction algorithm, so in practice it is worth choosing an optimal $M$ based on the experimental accuracy of phases control. We believe that extending the algorithm to consider the absolute values of the phases can help in determining the phase control errors. This is the subject of further research.

Finally, note that any interferometer architecture (for example, the Reck [12] or Clements [13] architecture) can be brought to the layered architecture considered here by adding controllable phase shifters in some channels. However, the possibility to directly apply our method to an arbitrary architecture remains unclear and is also the subject of further research.