Low-Depth Optical Neural Networks

TL;DR

This paper introduces a low-depth optical neural network architecture called OCTOPUS, which significantly improves noise robustness and scalability for machine learning tasks by compressing circuit depth logarithmically.

Contribution

The authors propose a novel low-depth ONN architecture that scales logarithmically with circuit depth, enhancing noise robustness and practical applicability.

Findings

LD-ONN shows exponential gain in noise robustness.

Numerical results on Letter Recognition dataset demonstrate improved stability.

OCTOPUS architecture can be used as a linear perceptron for classification.

Abstract

Optical neural network (ONN) is emerging as an attractive proposal for machine-learning applications, enabling high-speed computation with low-energy consumption. However, there are several challenges in applying ONN for industrial applications, including the realization of activation functions and maintaining stability. In particular, the stability of ONNs decrease with the circuit depth, limiting the scalability of the ONNs for practical uses. Here we demonstrate how to compress the circuit depth of ONN to scale only logarithmically, leading to an exponential gain in terms of noise robustness. Our low-depth (LD) ONN is based on an architecture, called Optical CompuTing Of dot-Product UnitS (OCTOPUS), which can also be applied individually as a linear perceptron for solving classification problems. Using the standard data set of Letter Recognition, we present numerical evidence showing that LD-ONN can exhibit a significant gain in noise robustness, compared with a previous ONN proposal based on singular-value decomposition [Nature Photonics 11, 441 (2017)].

Tables

Table 1. Table 1: Summary of the cost per layer for different ONN structures. All values corresponds to the order, O ​ ( ⋅ ) 𝑂 ⋅ O(\cdot) , for realizing one layer of the neural network with input dimension N 𝑁 N and output dimension M 𝑀 M .

	SVD	TLD	RLD
Circuit depth	$\max (N,M)$	$\log N$	$\log N$
Error scaling	$\max (N,M)$	$\log N$	$\log N$
Number of elements	$\max (N^2,M^2)$	$NM$	$N^{2}M$

Equations

$$\mathbf{y}=W\mathbf{x}.$$

$$\mathcal{D}(\mathbf{y},\tilde{\mathbf{y}})\equiv 1-\frac{\mathbf{y}\cdot \tilde{\mathbf{y}}}{\Vert \mathbf{y}\Vert \Vert \tilde{\mathbf{y}}\Vert },$$

$$\displaystyle f_{\bm{w}}(\bm{x})=\left\{\begin{array}[]{rcl}1&&\bm{w}\bm{\cdot x}\geqslant 0\\ 0&&\bm{w}\bm{\cdot x}<0\\ \end{array}\right.\ ,$$

$${\mathbf{v}}_1=\Gamma {\mathbf{v}}_{\mathbf{\alpha }}=[{\gamma }_1{\alpha }_1,{\gamma }_2{\alpha }_2,\cdots ,{\gamma }_N{\alpha }_N{]}^T.$$

$$\left[\begin{array}{c}{v}_{-}^{\mathrm{out}}\\ v_{+}^{\mathrm{out}}\end{array}\right]=H\left[\begin{array}{c}{v}_1^{\mathrm{in}}\\ v_2^{\mathrm{in}}\end{array}\right]=\frac{1}{\sqrt{2}}\left[\begin{array}{cc}1& -1\\ 1& 1\end{array}\right]\left[\begin{array}{c}{v}_1^{\mathrm{in}}\\ v_2^{\mathrm{in}}\end{array}\right].$$

$$T_i^{+}=\frac{1}{\sqrt{2}}\left[\begin{array}{cccccccc}1& 1& & & & & & \\ & & 1& 1& & & & \\ & & & & \ddots & & & \\ & & & & & & 1& 1\end{array}\right].$$

$$T_i^{-}=\frac{1}{\sqrt{2}}\left[\begin{array}{cccccccc}1& -1& & & & & & \\ & & 1& -1& & & & \\ & & & & \ddots & & & \\ & & & & & & 1& -1\end{array}\right].$$

$${\mathbf{v}}_{\mathbf{j}+1}=T_{n-j+1}^{+}{\mathbf{v}}_{\mathbf{j}}.$$

$${\mathbf{v}}_{\mathrm{out}}=T_1{T}_2\cdots T_n{\mathbf{v}}_1=\frac{1}{\sqrt{N}}\sum _{j=1}^N{\gamma }_j{\alpha }_j=\frac{1}{\sqrt{N}}\beta ,$$

$$W_j^{(i+1)}=\left[W_{j,1}^{(i+1)},W_{j,2}^{(i+1)},\cdots ,W_{j,N_i}^{(i+1)},W_{j,N_i+1}^{(i+1)}\right].$$

$$W^{(i)}={\left[W_1^{(i)},W_2^{(i)},\cdots ,W_{N_{i+1}}^{(i)}\right]}^T.$$

$${\mathbf{h}}^{(\mathbf{i}+1)}=F\left(W^{(i)}{\mathbf{h}}^{\prime (\mathbf{i})}\right),$$

$$S_j^n({\mathbf{v}}_1)=\{\begin{array}{rlrl}& T_{n-j+1}^{-}{\mathbf{v}}_{\mathbf{j}},& & j=1,\cdots ,n-1\\ & H{\mathbf{v}}_{\mathbf{n}},& & j=n.\end{array}$$

$$H\to \tilde{H}=\frac{1}{\sqrt{2}}\left[\begin{array}{cc}1+\delta & -1+\delta \\ 1-\delta & 1+\delta \end{array}\right].$$

$${\mathbf{v}}_{\mathrm{out}}=C\mathbf{y}=CW\mathbf{x},$$

$${\mathbf{v}}_{\mathrm{out}}^{\prime }=C{\mathbf{y}}^{\prime }=C\tilde{W}\mathbf{x}.$$

$$\mathbf{\delta y}\equiv [\delta y_1,\delta y_2,\cdots ,\delta y_M{]}^T\equiv ({\mathbf{y}}^{\prime }-\mathbf{y}).$$

$${\tilde{y}}_j$$

$$=y_i+\delta y_i^{\mathrm{re}}+O\left(\frac{(\delta y_i^{\mathrm{im}}{)}^2}{y_i}\right).$$

$${\tilde{y}}_j=y_j+\delta y_j,\delta y_j\in \mathbb{R}.$$

$$\mathcal{D}(\mathbf{y},\tilde{\mathbf{y}})$$

$$=1-\frac{\sum y_i{\tilde{y}}_i}{\sqrt{{\sum }_i{y}_i^2}\sqrt{{\sum }_i{\tilde{y}}_i^2}}$$

$$=1-\left(1-\frac{\sum y_i\delta y_i+\delta y_i^2/2)}{\Vert \mathbf{y}{\Vert }^2}\right)\left(1+\frac{\sum y_i\delta y_i}{\Vert \mathbf{y}{\Vert }^2}\right)+O\left(\frac{1}{\Vert \mathbf{y}{\Vert }^4}\right)$$

$$=\frac{1}{2\Vert \mathbf{y}{\Vert }^2}\sum \delta y_i^2+O\left(\frac{1}{\Vert \mathbf{y}{\Vert }^4}\right)$$

$$\sim \frac{N{\sigma }_y^2}{2N\overline{y_i^2}}$$

$$\sim \frac{{\sigma }_y^2}{\overline{y_i^2}}.$$

$$\Gamma \to \tilde{\Gamma }=\mathrm{diag}({\gamma }_1+\delta {\gamma }_1,{\gamma }_2+\delta {\gamma }_2,\cdots ,{\gamma }_N+\delta {\gamma }_N),$$

$$T_i\to {\tilde{T}}_i=\frac{1}{\sqrt{2}}\left[\begin{array}{cccccccc}1+{\delta }_1^{(i)}& 1-{\delta }_1^{(i)}& & & & & & \\ & & 1+{\delta }_2^{(i)}& 1-{\delta }_2^{(i)}& & & & \\ & & & & \ddots & & & \\ & & & & & & 1+{\delta }_{2^{i-1}}^{(i)}& 1-{\delta }_{2^{i-1}}^{(i)}\end{array}\right],$$

$$\tilde{T}={\tilde{T}}_1{\tilde{T}}_2\cdots {\tilde{T}}_n.$$

$$\tilde{T}=\frac{1}{\sqrt{N}}\left[\begin{array}{cccc}1+{\Delta }_1& 1+{\Delta }_2,& \cdots ,& 1+{\Delta }_N\end{array}\right]+O({\sigma }_I^2)$$

$${\Delta }_j=\sum _{i=1}^n{\delta }_j^{\prime (i)},$$

$$\displaystyle\delta^{\prime(i)}_{j}=\left\{\begin{array}[]{rcl}\delta^{(i)}_{\lfloor i/2^{n-j+1}\rfloor}&&\mod(i/2^{n-j},2)=0\\[6.0pt] -\delta^{(i)}_{\lfloor i/2^{n-j+1}\rfloor}&&\mod(i/2^{n-j},2)=1.\end{array}\right.$$

$${\tilde{\mathbf{v}}}_{\mathbf{n}}=\tilde{T}\tilde{\Gamma }{\mathbf{v}}_{\mathbf{\alpha }}\approx \frac{1}{\sqrt{N}}\left[\beta +\sum _{j=1}^N{\alpha }_j\delta {\gamma }_j+{\alpha }_j{\gamma }_j{\Delta }_j\right].$$

Full text

Introduction

Photonic computation represents an emerging technology enabling high-speed information processing with low energy consumption Flamini et al. (2018). Such technology can potentially be applied to solve many problems of machine learning, which has already created a significant impact on the physics community Biamonte et al. (2017); Mehta et al. (2018); Carleo and Troyer (2017); Carrasquilla and Melko (2017); Ma and Yung (2018); Bukov et al. (2018); Yang et al. (2018); Zhang et al. (2018); Gao et al. (2018). In particular, efforts have been made for decades in developing optical neural networks (ONNs) with different approaches Wagner and Psaltis (1987); Jutamulia and Yu (1996); Shen et al. (2017); Tait et al. (2017); Lin et al. (2018); Chang et al. (2018); Hughes et al. (2018); Bagherian et al. (2018); Penkovsky et al. (2019); Feldmann et al. (2019). Recently, much progress has been made in developing scalable on-chip photonic circuits Harris et al. (2017); Flamini et al. (2018); Wang et al. (2017); Carolan et al. (2015); Spring et al. (2013), leading to a new avenue towards large-scale implementation of ONNs. Compared with its free-space counterpart, on-chip ONN has advantages in terms of programmability and integrability Shen et al. (2017). This unconventional hardware architecture could potentially revolutionize the field of AI computing.

In order to achieve scalable ONNs, various circuit designs have been proposed recently Shen et al. (2017); Hughes et al. (2018); Bagherian et al. (2018), and they share similar characteristics, such as the scaling complexity of the circuit depth and the form of the multiport interferometers. In particular, ONN-based deep learning has been experimentally demonstrated Shen et al. (2017), by applying singular-value decomposition (SVD) for constructing any given linear transformation. Physically, these unitary transformations can be achieved with multiport interferometers Reck et al. (1994); Clements et al. (2016), together with a set of diagonal attenuators.

However, the circuit structure of SVD-ONN is only applicable for linear transformation represented by a square matrix; with $N$-dimensional input and $M$-dimensional output of data, the SVD approach of ONN requires $O(\max(N,M))$ layers of interferometers. As each layer will introduce errors to its output, the scalability of the SVD approach of ONN is limited by the errors scaling as $O(\max(N,M))$. Moreover, for machine-learning tasks of practical interest, both cases $N\gg M$ (e.g. image recognition He et al. (2016)) and $M\gg N$ (e.g. generative model Goodfellow et al. (2014)) are very common. Therefore, the SVD approach would require appending lots of ancillary modes to “square the matrix”, increasing the spatial complexity of the ONN.

To surmount the problem of robustness and flexibility, we propose an alternative approach of ONN for performing machine-learning tasks. Our ONN is constructed by connecting basic optical units, called Optical Computation of dot-Product Units (OCTOPUS), which optically outputs the dot-product of two vectors; the resulting circuit depth scales logarithmically $O(\log N)$. Even a single OCTOPUS can be applied as an optical linear perceptron Freund and Schapire (1999). In addition, the noise robustness of the OCTOPUS exhibits an exponential advantage compared with the SVD approach (see sm for the theoretical analysis ).

On the other hand, for constructing a deep neural network, we propose two variants of low depth ONN, called tree low depth (TLD) and recursive low depth (RLD) ONN. Both architectures involve OCTOPUS as basic optical computing units, and they are applicable for non-square transformation at each layer as well, i.e., $N\neq M$. The TLD-ONN requires fewer optical elements, but may cost more energy; the RLD-ONN involves a more complex structure, but it is more energy efficient. In terms of noise robustness, our numerical simulation suggests that TLD- and RTD-ONN have the same level of robustness, but both of them are significantly better than SVD-ONN.

Results

Linear Transformation— Given a one dimensional real vector $\bm{x}$ and an $N\times M$ real transformation matrix $W$, our goal is to optically achieve the following linear transformation

[TABLE]

In the SVD approach Shen et al. (2017), $N=M$ is assumed. Otherwise one needs to manually append many “[math]”s to square the corresponding matrix and vectors. Then, the matrix is decomposed as $W=U\Sigma V^{{\dagger}}$ (see Fig. 1), where $U$ and $V^{\dagger}$ are unitary matrices, and $\Sigma$ is a diagonal matrix. In optical implementation, $U$ and $V^{\dagger}$ can be realized with multiport interferometers of circuit depth $O(N)$, and $\Sigma$ can be realized with a set of attenuators or amplifiers Shen et al. (2017); Steinbrecher et al. (2018).

In contrast, our OCTOPUS solves the same problem by calculating the elements of the output $\bm{y}$ “one by one”. We require $M$ copies of the input $\bm{x}$, each of which serves as the input of one OCTOPUS. The $i$th OCTOPUS encodes the $i$th row of the matrix $W$ (denoted with $w_{i}$), and aims at calculating the $i$th element of the result $y_{i}=\bm{w_{i}}^{T}\cdot\bm{x}$.

The structure of OCTOPUS is shown in Fig. 2a (see Methods for details). At Part I, $\bm{w_{i}}^{T}$ is encoded with a set of tunable attenuators and phase shifters. In particular, the attenuators encode the magnitude of $\bm{w_{i}}^{T}$, while the phase shifters conditionally add a $\pi$ phase to the signal when the elements are negative. At Part II, we require a set of optical Hadamard transformation Reck et al. (1994); Clements et al. (2016). After each Hadamard transformation, we only trace the output port corresponding to the “sum” of its input (other paths denoted with dash lines are discarded). They are constructed as an interferometer tree of $n$ layers with totally $N=2^{n}$ input ports and $1$ output port. The amplitude of the final output becomes $\frac{1}{\sqrt{N}}\bm{w_{i}}^{T}\cdot\bm{x}=\frac{1}{\sqrt{N}}y_{i}$, which is the desired outcome multiplying a constant.

Note that the circuit depths of the SVD approach and the OCTOPUS approach are very different. The former scales linearly $O(N)$ and the latter logarithmically $O(\log N)$, leading to a dramatic difference in terms of the noise robustness against encoding errors of the optical elements. Specifically, let us denote the output vector subject to noise by $\bm{\widetilde{y}}$. The error can be quantified by cosine distance

[TABLE]

which has been widely adopted in classification problems Nair and Hinton (2010); Nguyen and Bai (2010); Dehak et al. (2011).

Our simulation results comparing the SVD approach with OCTOPUS against Gaussian noises are shown in Fig. 3a. For the SVD approach, the error $\mathcal{D}(\bm{y},\bm{\widetilde{y}})$ increases linearly with $N$. On the other hand, the error for OCTOPUS grows only very slowly. These results are consistent with the scaling of the circuit depths of the two approaches. In Supplemental Material sm , we provide further theoretical analysis on the noise effect. For SVD approach, the error scales linearly with the data size $N$, $\mathcal{D}(\bm{y},\bm{\widetilde{y}})\sim\sigma_{I}^{2}N+\sigma_{A}^{2}$, where $\sigma_{I}$ and $\sigma_{A}$ represent noise level for interferometers and attenuators respectively. However, for OCTOPUS, the error scales only logarithmically, $\mathcal{D}(\bm{y},\bm{\widetilde{y}})\sim\sigma_{I}^{2}\log{N}+\sigma_{A}$. This exponential advantages of the OCTOPUS approach agrees well with our numerical results (see Fig. 3a).

Linear perceptron— As OCTOPUS enables one to address each element of the input vector directly, it is possible to optically realize linear perceptrons Rosenblatt (1957); MacKay and Mac Kay (2003) with OCTOPUS for solving classification problems in machine learning.

For the binary case, the goal of linear perceptron is to output a hyperplane separating two classes of data labeled by either [math] or $1$, allowing us to make prediction on the unlabelled new data. More precisely, with a set of training data $\{\bm{x}\}$, one needs to determine the parameters $\bm{w}$ for the following function $f_{\bm{w}}(\bm{x})$:

[TABLE]

which can be realized with a single use of OCTOPUS followed by an appropriate measurement. As shown in Fig. 2d, at the $i$th iteration, we use training data $\bm{x^{(i)}}$ as the input, determining whether its corresponding label $y^{(i)}=f_{\bm{w}}(\bm{x^{(i)}})$. If it is not, the weight is updated according to $\bm{w}\leftarrow\bm{w}+\eta\bm{x^{(i)}}$, where $\eta$ is the learning rate.

Fig. 3b shows the simulation results of our optical linear perceptron with “Iris Data Set” Blake (1998). We define error rate as the rate of providing incorrect prediction on the label. Remarkably, the error rate remains under $0.1$ when the noise level $\sigma\leqslant 0.05$.

Low-depth ONN (LD-ONN)— It is well known that linear perceptron performs well with relatively simple tasks. However, for problems involving complicated non-linear relations, one may consider deep neural networks. In the following, we present two variants of multi-layered ONN, namely, Tree Low-Depth (TLD) and Recursive Low-Depth (RLD) ONN. Both of them share a similar structure, as illustrated in Fig. 4a. Again, the input data $\bm{x}$ is encoded at the first layer containing an attenuator and a phase shifter at each node. Optical computation is performed at each neuron (denoted with colored circle), encapsulating the trainable parameters of the networks. Furthermore, the optical paths at the top of each layer represent the “biases” of the corresponding layer.

As shown in Fig. 4b,c, for TLD-ONN, each neuron consists of an OCTOPUS together with a nonlinear activation function, which can be physically realized with non-linear crystal Shen et al. (2017), measurement Hughes et al. (2018) or optical amplifier Connelly (2007); sm . Then, each path is distributed uniformly to many paths, which are the inputs of the neurons at the next layer. More details are given in the Methods section and Supplemental Material sm .

Note that in TLD-ONN, the OCTOPUS only picks one path as its output; many paths are discarded. To realize a deep ONN, one may need a strong light source or amplify the signal at each layer. Alternatively, we may re-structure the ONN, which is the motivation for developing the RLD-ONN approach.

As shown in Fig. 4b and c, for RLD-ONN, the input signals first pass through a set of trainable attenuators and phase shifters, followed by a ($3$ or $2$-layer) interferometer tree. Different from TLD-ONN, no signal are discarded after these steps. At this point, only two paths (such as “$1$” and “$2$” in Fig. 4c) are fully connected to the corresponding input ports. In order to connect all other output paths, the remaining paths are sent to interferometer trees with smaller size recursively. Then, the nonlinear activation is applied to all output paths. If the number of neurons at the next layer is larger than the current number of output paths, the output paths can be expanded with $50/50$ beam splitter (Fig. 4b); if it is less, one can just discard several output paths (Fig. 4c)not . Note that the way of connecting input and output paths are not unique, so further optimization can be performed. Although RLD-ONN requires more optical elements, the circuit depth remains to be logarithmic.

The universality approximation theorem states that a feed-forward neural network with as few as a one hidden layer can approximate any continuous function to an arbitrary accuracy (under some mild assumptions on the activation function) Hornik (1991), which is the foundation of neural computation. In Methods section, we show that the transformation of TLD-ONN is equivalent to standard feed-forward nerual network. And in sm we show that for any given from of one-hidden-layer TLD-ONN, there always exist a RLD-ONN that is equivalent to it. An illustration is also provided in Fig. S2. Therefore, both TLD- and RLD-ONN proposed in this work are universal.

To compare the performance of LD-ONNs with SVD-ONN Shen et al. (2017), we perform numerical simulation on the “Letter Recognition” data set Blake (1998), classifying letters “A”, “B”, “C” and “D” (see Methods section for technical details). The ONN used in the simulation contains one hidden layer with $64$ neurons. While the SVD- and TLD-ONN are trained with standard back-propagation method, the RLD-ONN is trained with “forward propagation”Shen et al. (2017). In this work, we consider the training as a pre-processing, i.e., the parameters of the network are first trained in conventional computer. However, the training can also be realized optically with little assistance from electronic devices sm . As shown in Fig. 4d, when noise level $\sigma=0$, the error rates for SVD-, TLD- and RLD-ONN are $6.4\%$, $6.4\%$, $6.6\%$ respectively. These values may be improved by further optimizing the hyperparameters. As the noise level increases, TLD- and RLD-ONN have comparable error rates, but both of them are significantly lower than the error rate of SVD-ONN.

Discussion

A summary of SVD-, TLD- and RLD-ONN is given in Table. 1, providing a comparison of the cost for an ONN layer with input dimension $N$ and output dimension $M$. Both the TLD- and RLD-ONN have logarithmic circuit depth, leading to exponential improvements on the error scalings compared with the SVD approach. Furthermore, TLD-ONN requires less number of optical elements, but we note that it also requires discarding more paths during the implementation. On the other hand, RLD-ONN requires discarding much fewer paths (same as SVD-ONN), but at the cost of a larger number of optical elements.

Note that several simplifications can further be made on the LD-ONN structures. Firstly, as discussed in Ref sm , the ONN can be binarized with amplifier working at the saturation regime. With the binarization of the weight, biases and activation function, the attenuators (for magnitude encoding) at each OCTOPUS can be removed. Secondly, instead of encoding the parameters at the phase shifter and attenuators, they can also be encoded at the interferometers. More specifically, it is possible to remove the part I of OCTOPUS, and replace the Hadamard transformation at Part II by tunable interferometers. Thirdly, after the training process, the network can be “compressed” with the “pruning” techniqueHan et al. (2015): removing all paths with weights below a threshold. Above improvements or revision could helps reducing the complexity of the hard-ware architectures.

To conclude, based on OCTOPUS, we present a new architecture of ONN for machine learning, which provides exponential improvements on the robustness against encoding error. We discussed different schemes of optical linear transformation, linear perceptron, and two variants of multi-layered ONNs. Numerical simulations with random transformations and standard machine learning data sets are employed to justify the robustness of our schemes. The proposed LD-ONN can be directly implemented with current photonic circuitsFlamini et al. (2018) technology. Our proposal, combined with appropriate realization of nonlinear activation (For example, the scheme in Ref Williamson et al. (2019)), provides a possible solution to solving machine-learning tasks of industrial interest with robust, scalable, and flexible ONNs.

Methods

OCTOPUS. We consider two vectors $\bm{\alpha}=\left[\alpha_{1},\alpha_{2},\cdots,\alpha_{N}\right]^{T}$ and $\bm{\gamma}=\left[\gamma,\gamma_{2},\cdots,\gamma_{N}\right]^{T}$ with $\alpha_{i}\in\mathbb{R}$ and $\gamma_{i}\in[-1,1]$, and assume $N=2^{n}$ with $n\in\mathbb{Z}^{+}$. They corresponds to $\bm{x}$ and $\bm{w_{i}^{T}}$ in Fig. 2a respectively. The OCTOPUS aims at calculating $\beta=\bm{\gamma\cdot\alpha}$. The input signal has $N$ paths, and the amplitude at the $i$th path is set to be $\alpha_{i}$, so it can simply be denoted with $\bm{v_{\alpha}}=\bm{\alpha}$.

At part I, $\bm{\gamma}$ is encoded with a set of attenuators and phase shifters at each path. At the $i$th path, the attenuator controls the magnitude of $\gamma_{i}$, while the phase shifter determines the sign of $\gamma_{i}$ ( when $\gamma_{i}$ is negative, the phase shifter performs a $\pi$ shift on the input signal). The total transformation can be represented by the diagonal matrix $\Gamma={\rm diag}(\gamma_{1},\gamma_{2},\cdots,\gamma_{N})$. So the output of part I can be represented by

[TABLE]

At Part II, all elements in Eq. (6) are “summed up” with an $n$-layer interferometer tree. The $j$th layer of the tree contains $2^{n-j}$ interferometers, each of which performs the Hadamard transformation on two nearest neighbour paths. If we denote the input of an interferometer as $[v^{\rm in}_{1},v^{\rm in}_{2}]^{T}$, and the output as $[v^{\rm out}_{-},v^{\rm out}_{+}]^{T}$, each interferometer performs the Hadamard transformation Clements et al. (2016) at two nearest neighbour paths

[TABLE]

We denote the input of the $j$th layer of the tree as $\bm{v_{j}}$, and separate the output paths at this layer into two groups. The “$+$” group contains all output paths corresponding to $v^{+}_{\rm out}$, while the “$-$” group contains those corresponding to $v^{-}_{\rm out}$. They can be represented by $T_{n-j+1}^{+}\bm{v_{j}}$ and $T_{n-j+1}^{-}\bm{v_{j}}$ respectively, where $T_{i}^{+}$ and $T_{i}^{-}$ are the following $2^{i}\times 2^{i-1}$ matrices:

[TABLE]

[TABLE]

[TABLE]

To obtain the summation, we discard the “$-$” group, and take “$+$” group as the input of the next [$(j+1)$th] layer

[TABLE]

So the total transformation at part II can be represented by:

[TABLE]

which is just the dot-product result $\beta$ multiplying a constant $\frac{1}{\sqrt{N}}$.

Amplitude measurement. To obtain the optical computation results, one should extract both the magnitude and the sign of the output. Here, we provide a measurement scheme in Fig. 2b. Suppose the amplitude of the signal to be measured is $\beta$, firstly, the signal is splitted into two paths with amplitude $\beta_{\rm m}$ for main path and $\beta_{\rm anc}$ for ancillary path, where $|\beta_{\rm anc}|\ll|\beta_{\rm m}|$ is required. By measuring the intensity at the main path, one obtains the magnitude of $\beta$. Then, we introduce a reference path with signal amplitude $\beta_{\rm ref}$, and interfere it with the ancillary path. We require that $\beta_{\rm ref}\gtrsim\beta_{\rm anc}$. The ancillary path and the reference path are interfered with a $50/50$ beam splitter. Obviously, if the intensity is enhanced after the interference, the sign of $\beta$ should be “$+$”, otherwise the sign should be “$-$”.

TLD-ONN structure. As can be seen in Fig. 4, if there are $N_{i}$ neurons and $N_{i+1}$ neurons at the $i$th and $(i+1)$th layer, there are totally $N_{i}\times N_{i+1}$ output paths from the $i$th layer. Actually, the signal amplitudes of output paths coming from the same neuron (totally $N_{i+1}$ paths) are identical. So we represent the output signals from the $i$th layer with an $N_{i}$ dimensional vector $\bm{h^{(i)}}=\left[h^{(i)}_{1},h^{(i)}_{2},\cdots,h^{(i)}_{N_{i}}\right]^{T}$, where $h^{(i)}_{j}$ represents the amplitude of all paths output from the $j$th neuron of the $i$th layer. Note that the input of the network is just $\bm{h^{(0)}}$.

To introduce the bias at each layer, the input of the $(i+1)$th layer is the output of the $i$th layer adding an ancillary path (at the top of each layer in Fig. 4a) with constant amplitude. Without loss of generality, we assume this constant to be 1. Therefore, there are $(N_{i}+1)$ input paths for each neuron at the $(i+1)$th layer, which we represent with $\bm{h^{\prime(i)}}=\left[h^{(i)}_{1},h^{(i)}_{2},\cdots,h^{(i)}_{N_{i}},1\right]^{T}$.

Each neuron contains an OCTOPUS encoding a set of parameters at its part I. We denote all parameters corresponding to the $j$th neuron of the $(i+1)$th layer as

[TABLE]

Here, $W^{(i+1)}_{j,k}$ denotes the parameter encoded at the $(i+1)$th layer, the $k$th path of the $j$th neuron, and $W^{(i+1)}_{j,N_{i}+1}$ is the bias of the $(i+1)$th layer.

As illustrated in Fig.4b, c for TLD-ONN, at each neuron of the $(i+1)$th layer, the signal is first pass through OCTOPUS, and then a nonlinear activation $f(x)$ (see Supplemental Materials and Fig. S1), and is finally distributed uniformly to $N_{i+2}$ paths. We define the linear transformation matrix from the $i$th layer to the $(i+1)$th layer as

[TABLE]

With a little thought, one can find that the relation between $\bm{h^{(i+1)}}$ and $\bm{h^{(i)}}$ are given by

[TABLE]

where $F(x)=\frac{1}{\sqrt{N_{i+2}}}f\left(\frac{1}{\sqrt{N_{i}}}x\right)$ is the rescaled nonlinear activation function. As can be seen, the transformation at TLD-ONN is equivalent to the standard feed-forward neural network. The training process is discussed in the Supplemental Material sm .

RLD-ONN structure. The general structure of RLD-ONN is similar to LD-ONN (Fig. 4a). The difference lies in the transformation performed by each neuron. As shown in Fig. 4b, c, each neuron contains three parts: encoding, recursively connecting, and nonlinear activation. The encoding part is identical to part I of OCTOPUS, and the nonlinear activation is discussed in Supplemental Material sm . So we focus on the recursively connecting part in the middle. Generally, the goal of this part is to make all input connect to all output ports while maintaining the logarithmic circuit depth, and does not discard any paths. It turns out that this can be realized with a set of revised interferometer tree.

We recall from the Part II of OCTOPUS that at the $j$th layer of the interferometer tree, its input paths are represented with $\bm{v_{j}}$, and the output paths are separated into “$+$” group and “$-$” group. As shown in Fig. 4b, c, the “$+$” group, as usual, serves as the input of the next layer. But instead of discarding the “$-$” groups, they are also traced. For an $n$-layer interferometer tree with input signal $\bm{v_{1}}$, we denote its output at the $j$th layer with $S^{n}_{j}(\bm{v_{1}})$, which is given by

[TABLE]

For $j<n$, the output is just the “$-$” group at the corresponding layer; and for $j=n$, the output consists both “$+$” and “$-$” groups.

Obviously, only the output ports at the $n$ layer $S^{n}_{n}(\bm{v_{1}})$ are connected to all input ports. Signals from the remaining output ports only contain local information of the input they are connected to, so the network is not expected to work well if one use them directly. The key idea of the “recursive” is that the interferometer trees described by Eq. (14) are recursively appended until all input and output ports are fully connected. Take Fig. 4c as an example. We first apply a $3$-layer interferometer tree, after which the first and second paths are connected to all input ports. We then apply a $1$-layer tree to the 3rd and 4th paths, and $2$-layer tree to 5th-8th paths, after which only the A and B paths do not connect to all input ports. So we finally apply $1$-layer tree to A and B paths. Then, all output and input ports are fully connected. Formally, for input vector $\bm{x}$ of dimension $N=2^{n}$, the transformation performed by the recursive structure $\bm{y}=U(\bm{x})$ is described by Alg. 1. Moreover, the training method for RLD-ONN is discussed in details in sm .

Technical details for simulations. In the simulation, the errors are introduced to attenuators and interferometers unless specified. For attenuator encoding the value $w$, the error is introduced by making the replacement $w\rightarrow w+\delta w$; for each interferometer, the replacement is

[TABLE]

In general, the error terms should be complex values. But in Supplemental Material sm , we show that when $|\delta w|\ll w$ and $|\delta|\ll 1$, the effect of the imaginary parts of $\delta w$ and $\delta$ are an order smaller than the real part. So the noise effect can be well approximated by restricting $\delta w$ and $\delta$ to be real. We further assume that $\delta$ and $\delta w$ are Gaussian noises drawn from $\mathcal{N}(0,\sigma_{I}^{2})$ and $\mathcal{N}(0,w^{2}\sigma_{A}^{2})$ respectively, where $\sigma_{I}$ and $\sigma_{A}$ are noise level for interferometers and attenuators respectively.

In the simulation of linear perceptron (Fig. 3b), the input data of “Iris” data set are four-dimensional vectors describing different properties of a particular “Iris” flower. We test our algorithm on two kinds of Iris, “Setos” and “Versicolour”, and label them with “1” and “0” respectively. Totally 100 pairs of data and labels are used, which are separated to training set (with size 60) and testing set (with size 40).

In the simulation of deep ONN (Fig. 4d), the networks contain one hidden layer with 64 neurons, and we assume that the activation function is inverse square root unit (ISRU), which is realized by optical amplifier (see Supplemental Material sm ). For a given Letter, there are totally 16 primitive numerical attributes, so the input data is $16$ dimensional. Each elements of the input data are rescaled to the interval $[0,1]$. The dimension of output vectors are four, and the desired output for each letter is $\bm{y_{a}}=[1,0,0,0]^{T}$ for “A”, $\bm{y_{b}}=[0,1,0,0]^{T}$ for “B”, $\bm{y_{c}}=[0,0,1,0]^{T}$ for “C” and $\bm{y_{d}}=[0,0,0,1]^{T}$ for “D” respectively. We use in total 2880 pairs of inputs and labels, which are separated to training set (with size 1920) and testing set (with size 960).

The training process follows Algorithm.S1 and Algorithm.S2 in sm . For SVD approach, after obtaining the well-trained $W^{(i)}$, we append it with enough “[math]” elements until they become square matrices, and the hidden layers, input, output vectors are appended with “[math]” accordingly. Singular-Value decomposition is performed at the matrices $W^{(i)}$, after which the parameters for attenuators and the M-Z interferometers are be determined with the method given in Ref Clements et al. (2016). To make the comparison fair enough, we set the learning rate ($0.01$), mini-batch size ($32$) and initial guess of all trainable parameter (drawn from $[-0.1,0.1]$ with uniform probability) to be the same for all approaches. The results shown in Fig. 4d are averaged over 10 runs with different generated random noise.

I Error estimation for linear transformation

Suppose we are given an input vector $\bm{x}$ and the transformation matrix $W$, our goal is to comput $\bm{y}=W\bm{x}$ optically. We restrict that $\bm{x}\equiv[x_{1},x_{2},\cdots,x_{N}]^{T}$, $\bm{y}\equiv[y_{1},y_{2},\cdots,y_{M}]^{T}$ and the transformation matrix $W$ to be real, as it is the common scenario of machine learning applications. Ideally, the amplitude of output signal of either OCTOPUS or SVD approaches are given by $\bm{v_{\rm out}}$, which satisfies

[TABLE]

where $CW$ is the transformation performed by the photonic circuits. As discussed in the Method section of [Smai, ], for OCTOPUS approach, the constant $C=1/\sqrt{N}$; for SVD, the constant $C=1$. When the encoding error is introduced, the imperfect transformation is replaced by $\widetilde{W}$, and the above equation becomes

[TABLE]

We define the error term as

[TABLE]

It is in general a complex vector, i.e., $\delta y_{j}=\delta y_{j}^{\rm re}+i\delta y_{j}^{\rm im}$ with $\delta y_{j}^{\rm re},\delta y_{j}^{\rm im}$ to be real and nonzero. While the vector $\bm{v^{\prime}_{\rm out}}$ is the amplitude of output signal, the computation result, $\bm{\widetilde{y}}\equiv[\widetilde{y}_{1},\widetilde{y}_{2},\cdots,\widetilde{y}_{M}]^{T}$, can only be estimated after the measurement (see Methods in [Smai, ]). If we assume the measurement process to be ideal, each element is given by

[TABLE]

Even though the imaginary part of the error is in general comparable to the real part in real experimental implementation, when $|\delta y_{i}|\ll y$, the effect of imaginary part to the final result is negligible compared with the real part. Therefore, when estimating the noise effect, we can safely simplify the noise model as

[TABLE]

For similar reason, in the following discussion, all encoding errors are assumed to be real.

I.1 Cosine Distance

Without loss of generality, we can assume the error term $\bm{\delta y}$ has zero mean and variance $\sigma_{y}^{2}\equiv\text{Var}(\delta y_{i})=\text{Var}(\tilde{y}_{i})$. We use cosine distance between $\bm{y}$ and $\bm{\widetilde{y}}$ to quantify the effect of error. When $\bm{y}\gg\bm{\delta y}$, it can be estimated by

[TABLE]

In the following, we will estimate the value of $\sigma_{y}^{2}$ for OCTOPUS and SVD approaches separately.

I.2 Error for OCTOPUS approach

We should first estimate the error of a single OCTOPUS. We denote the transformation performed by imperfect OCTOPUS with encoding error as $\widetilde{\beta}=\bm{\widetilde{\gamma}\cdot\alpha}$.

We introduce encoding error for attenuators at part I:

[TABLE]

and interferometers at part II:

[TABLE]

where $\delta^{(i)}_{k}$ and $\delta\gamma_{j}$ and are assumed to be Gaussian noise drawn from $\mathcal{N}(0,\sigma_{I}^{2})$ and $\mathcal{N}(0,\gamma_{j}^{2}\sigma_{A}^{2})$ independently. The total transformation at Part II is replaced by

[TABLE]

After some calculation, we obtain

[TABLE]

where

[TABLE]

and

[TABLE]

So the output of part II becomes:

[TABLE]

Therefore, the estimated value of $\beta$ is

[TABLE]

Since $\delta\gamma_{j}$ and $\delta^{\prime(i)}_{j}$ are independent of each other, for large enough $N$ we have

[TABLE]

and

[TABLE]

Therefore, the variant of $\widetilde{\beta}$ is estimated as

[TABLE]

Since each element $y_{i}$ are calculated with its corresponding OCTOPUS, according to Eq. (S-19), the variant of $y_{i}$ satisfies

[TABLE]

So $\mathcal{D}(\bm{y},\bm{\widetilde{y}})$ can be estimated as

[TABLE]

The error contributed from Part I does not increase as $N$ increase, while the error contributed from part II increase as $O(n)=O(\log N)$. These results are consistent with the circuit depth of both part I and part II, as well as the numerical results in Fig.3a of [Smai, ].

I.3 Error for SVD approach

To begin with, we first review how SVD approach realize arbitrary real value linear transformation with photonic circuit of totally $2N+1$ layers. We can denote the amplitudes at the $j$th layer and $i$th path as $x^{(j)}_{i}$, and use $\bm{x^{(j)}}=[x^{(j)}_{1},x^{(j)}_{2},\cdots,x^{(j)}_{N}]^{T}$ to represent the signal at the $j$th layer. The input ports and output ports of the ONN correspond to $\bm{x^{(0)}}=\bm{x}$ and $\bm{x^{(2N+1)}}=\bm{y}$ respectively.

We use $M^{(j)}$ to represent the transformation at the $j$th layer, so the transformation from $\bm{x^{(j)}}$ to $\bm{x^{(j+1)}}$ can be denoted as

[TABLE]

While $M^{(N+1)}$ is a diagonal matrix, $M^{(j\neq N+1)}$ are unitary matrices containing only the interaction between nearest neighbour paths. The interaction between nearest neighbor paths are introduced with M-Z interferometers. Since we restrict all transformation to be real, they can be represented as

[TABLE]

We denotes $R_{i}(\theta)$ as the transformation performed at the $i$th and $(i+1)$th paths while all other paths remain unchanged. If we assume $N$ is even, the transformation at $j$th layer can be represented by

[TABLE]

Here, $R_{k}(\theta_{2k-1,j})$ and $R_{k}(\theta_{2k,j})$ corresponds to the M-Z interferometer at the $j$th layer and connects the $k$th and $(k+1)$th paths; $s_{j}\in[-1,1]$ corresponds to the attenuator at the $(N+1)$th layer and the $j$th paths.

The total transformation is given by

[TABLE]

It has been shown that arbitrary real $W$ can be realized by choosing $\theta_{i,j}$ and $s_{i}$ appropriately [SClements et al., 2016].

To study the effect of encoding error, we can also do the following replacements

[TABLE]

Similarly, we assume that $\delta\theta_{i,j}$ and $\delta s_{i}$ are real, and are drawn from normal distribution $\mathcal{N}(0,\sigma_{I}^{2})$ and $\mathcal{N}(0,s_{i}^{2}\sigma_{A}^{2})$ respectively. The transformation at each layer then becomes:

[TABLE]

where

[TABLE]

The final output then becomes

[TABLE]

where

[TABLE]

Since $M^{(i)}$ are either unitary matrix or diagonal matrix with the values restricted to $[-1,1]$, the order of $\bm{\delta x^{(j)}}$ is determined by $\delta M^{(j)}$. We recall that the variance of $\delta\theta_{i,j}$ and $\delta s_{i}$ are $\sigma_{I}^{2}$ and $s_{i}^{2}\sigma_{A}^{2}$ respectively. Let $\delta x^{(j)}_{i}$ to be the $i$th element of $\bm{\delta x^{(j)}}$, we have

[TABLE]

Since $\delta y_{i}=\sum_{j=1}^{2N+1}\delta x^{(j)}_{i}$ and $\delta x^{(j)}_{i}$ are independent of each other, for large $N$ we have

[TABLE]

Combining with Eq. (S-6), the cosine distance can be estimated as

[TABLE]

So $\mathcal{D}(\bm{y},\bm{\widetilde{y}})$ increases as $O(N)$, which agrees well to the linear depth of SVD circuit and the numerical result in Fig.3a of [Smai, ].

II Amplifier as activation function

There are many ways to generate different types of nonlinear activation, such as with the saturable absorbers [SShen et al., 2017] or via measurement [SHughes et al., 2018]. With the former one, the transformation is realized with the speed of light, but the signal will attenuate when passing through each layer. With the measurement approach, the intensity will not decrease as the network size increase, but the computation speed would be reduced by the measurement process. Here, we introduce an alternative approach based on optical amplifier, which does not suffer from signal strength reduction, and maintain the high computation speed. But the most appropriate solution to the activation function depends on the practical scenarios, and require further studies.

For an optical amplifier, the power of input and output signals are given by $P_{\rm in}=|x|^{2}$ and $P_{\rm out}=|f(x)|^{2}$ respectively, where $x$ and $f(x)$ are the amplitude of input and output signal. The gain of an optical amplifier, $g=P_{\rm out}/P_{\rm in}$, generally satisfies [SConnelly, 2007]:

[TABLE]

where $P_{s}$ is the saturation power. So for an amplifier that maintains the signal phase, the nonlinear relation between input and output is

[TABLE]

which is exactly the inverse square root unit (ISRU) [SCarlile et al., 2017]. Through out our simulation, Eq. (S-35) is used as the activation function. On the other hand, if the amplifier works in the saturation regime, i.e. $x\ll\sqrt{P_{s}}$, we have the following binarized activation function

[TABLE]

with which one can construct a binarized neural network [SCourbariaux et al., 2016; Hubara et al., 2016]. This type of neural networks only need the amplitudes $\bm{h^{(i)}}$ and the elements of matrices $W^{(i)}$ to be either $+1$ or $-1$, but still have comparable performance to traditional networks. So one may encode the parameters only without tunnable attenuator, and dramatically simplify the structures.

In Fig. S1, we show the input-output relations for different values of $P_{s}$, corresponding to Eq. (S-35) and Eq. (S-36) respectively.

III Training for TLD-ONN

The discrepancy between desired output (or label), $\bm{y}$ and the output of neural network $\bm{h^{(L)}}$ is quantified by the loss function $\mathcal{L}$, which is taken to be the mean square error in this work

[TABLE]

During training, our goal is to minimize the loss function for the given training set by tuning $W^{(i)}$ containing the weights and biases of the $i$th layer. The process follows the standard back-propagation method as shown in Algorithm. S1, where we randomly generate a minibatch containing $N_{b}$ pairs of input data and labels $(\bm{x_{i_{b}}},y_{i_{b}})\in\mathcal{B}$ at each iteration. The training can be realized solely in the electronic devices, after which the well-trained parameters are encoded to the ONN setup. Alternatively, one may train the network with the assistant of OCTOPUS: by executing all linear transformation steps (commented steps) with OCTOPUS, the training can be accelerated dramatically.

We also note that it is possible to implement the “forward propagation” approach, which train the ONN directly by tuning the attenuators and phase shifters, and obtain the gradient by directly measure the output of the ONN [SShen et al., 2017].

IV Training for RLD-ONN

Since RLD-ONN has a special structure, standard back-propagation is no-longer available. Instead, one can train the RLD-ONN with “forward-propagation” method, i.e., perturb each parameter directly, and update the network according to the gradient of the cost function with respected to all trainable parameters [SShen et al., 2017].

Similar to TLD-ONN, we denote $w^{(i)}_{j,k}$ as the trainable parameter (encoded at the attenuators and phase shifters) at the $i$th layer, the $j$th neuron, and the $k$th path. We simply use $\Theta$ to represent all parameters of the neural network, and use $G(\bm{x},\Theta)$ to represent the output of the neural network with input $\bm{x}$ and parameter $\Theta$. $G(\bm{x},\Theta)$ can be calculated in an electronic device, or it can be estimated directly at RLD-ONN. The loss function [Eq. (S-37)] then becomes

[TABLE]

We further define the perturbed output the neural network $G\left(\bm{x},\Theta,w^{(i)}_{j,k},\delta\right)$, which is the neural network output with parameter $\Theta$, except for the element $w^{(i)}_{j,k}$ changed as $w^{(i)}_{j,k}\leftarrow w^{(i)}_{j,k}+\delta$. The derivative of the loss function with respected to $w^{(i)}_{j,k}$ can be estimated with

[TABLE]

where we assume $\delta\ll w^{(i)}_{j,k}$. $\Theta$ is updated according to the gradient of the loss function. At each iteration $i_{b}$, we calculate the gradient with respected to a minibatch containing $N_{b}$ pairs of input data and label $(\bm{x_{i_{b}}},y_{i_{b}})\in\mathcal{B}$, and perform the gradient descendent base on it. The full training algorithm is given in Algorithm. S2.

V Universality of RLD-ONN

To show the universality of RLD-ONN, one just need to ensure that for TLD-ONNs with any given parameters, there are a RLD-ONN that can reduce to it. This turns out to be true. We consider a one-hidden-layer TLD-ONN with $N_{\text{i}}$ input, $N_{\text{h}}$ hidden, and $N_{\text{o}}$ output neurons respectively. Suppose the parameters encoded at the $i$th neuron of the hidden layer is

[TABLE]

and the parameters encoded at the $i$th neuron of the output layer is

[TABLE]

In the following, we show that this TLD-ONN is identical to an RLD-ONN with $N_{\text{i}}$ input, $N_{\text{h}}\times N_{\text{o}}$ hidden, and $N_{\text{o}}$ output neurons. We denote the parameters encoded at the $i$th neuron at the hidden layer and output layer as

[TABLE]

and

[TABLE]

Firstly, we let the first $N_{\text{h}}$ hidden neurons just connect to the first output layer, the $(N_{\text{h}}+1)$ to $2N_{\text{h}}$ hidden neurons just connect tot he second output and so on. In other words, at the output layer, we set

[TABLE]

As can be seen, there are totally $N$ nonzeros input at the output layers, which is the same as its corresponding TLD-ONN. Then, one just need to ensure that the input signals of the output layer are identical to those in TLD-ONN, by setting $\bm{V^{(h)}_{i}}$ appropriately. For the $i$th neuron at the hidden layer, the signal before entering the nonlinear activation can be represented by

[TABLE]

with $\bm{x^{\prime}}=[\bm{x}^{T},1]^{T}$ the input of ONNs. As discussed in the main text, $U$ is realized by a set of interferometers, so it corresponds to a unitary transformation. In order words, $U(\bm{x})=U\bm{x}$ for certain unitary matrix $U$. So we have

[TABLE]

for certain $\bm{u_{j}}=[u_{j,1},u_{j,2},\cdots,u_{j,N}]^{T}$. As can be inferred from Eq. (S-44), only the path corresponds to $h_{i,\lceil i/N_{\text{o}}\rceil}$ are used in the next layer. Therefore, to mach the corresponding TLD-ONN, the only constrain is

[TABLE]

which can be satisfied by setting

[TABLE]

where we have defined $\bm{u_{j}^{-1}}=[u_{j,1}^{-1},u_{j,2}^{-1},\cdots,u_{j,N}^{-1}]^{T}$. An illustration of the case $N_{\text{i}}=2$,$N_{\text{h}}=3$,$N_{\text{i}}=2$ is shown in Fig. S2.

Here, we have shown that RLD-ONN with more hidden neurons can reduce to TLD-ONN, so it is universal. But in practical application, the RLD-ONN may need much less hidden neurons than $N_{\text{h}}\times N_{\text{o}}$ to achieve a comparable performance with TLD-ONN.