Noise-Resilient Designs for Optical Neural Networks
All analog signal processing is fundamentally subject to noise, and this is
also the case in modern implementations of Optical Neural Networks (ONNs).
Therefore, to mitigate noise in ONNs, we propose two designs that are
constructed from a given, possibly trained, Neural Network (NN) that one wishes
to implement. Both designs have the capability that the resulting ONNs gives
outputs close to the desired NN.
  To establish the latter, we analyze the designs mathematically. Specifically,
we investigate a probabilistic framework for the first design that establishes
that the design is correct, i.e., for any feed-forward NN with Lipschitz
continuous activation functions, an ONN can be constructed that produces output
arbitrarily close to the original. ONNs constructed with the first design thus
also inherit the universal approximation property of NNs. For the second
design, we restrict the analysis to NNs with linear activation functions and
characterize the ONNs' output distribution using exact formulas.
  Finally, we report on numerical experiments with LeNet ONNs that give insight
into the number of components required in these designs for certain accuracy
gains. We specifically study the effect of noise as a function of the depth of
an ONN. The results indicate that in practice, adding just a few components in
the manner of the first or the second design can already be expected to
increase the accuracy of ONNs considerably.
Gianluca Kosmella, Ripalta Stabile, Jaron Sanders
# Noise-Resilient Designs for Optical Neural Networks

###### Abstract

All analog signal processing is fundamentally subject to noise, and this is also the case in modern implementations of Optical Neural Networks (ONNs). Therefore, to mitigate noise in ONNs, we propose two designs that are constructed from a given, possibly trained, Neural Network (NN) that one wishes to implement. Both designs have the capability that the resulting ONNs gives outputs close to the desired NN.

To establish the latter, we analyze the designs mathematically. Specifically, we investigate a probabilistic framework for the first design that establishes that the design is correct, i.e., for any feed-forward NN with Lipschitz continuous activation functions, an ONN can be constructed that produces output arbitrarily close to the original. ONNs constructed with the first design thus also inherit the universal approximation property of NNs. For the second design, we restrict the analysis to NNs with linear activation functions and characterize the ONNs' output distribution using exact formulas.

Finally, we report on numerical experiments with LeNet ONNs that give insight into the number of components required in these designs for certain accuracy gains. We specifically study the effect of noise as a function of the depth of an ONN. The results indicate that in practice, adding just a few components in the manner of the first or the second design can already be expected to increase the accuracy of ONNs considerably.

keywords: Optical Neural Networks, Law of Large Numbers, Universal Approximation +
Footnote †: journal: Journal of LaTeX Templates

## 1 Introduction

Machine Learning (ML) is a computing paradigm in which problems that are traditionally challenging for programmers to explicitly write algorithms for, are solved by learning algorithms that improve automatically through experience. That is, they "learn" structure in data. Prominent examples include image recognition [1], semantic segmentation [2], human-level control in video games [3], visual tracking [4], and language translation [5].

Classical computers are designed and best suited for serialized operations (they have a central processing unit and separated memory), while the data-driven ML approach requires decentralized and parallel calculations at high bandwidth as well as continuous processing of parallel data. To illustrate how ML can benefit from a different architecture, we can consider performance relative to the number of executed operations, also indicated as Multiply--Accumulate Operation (MAC) rates, and the energy efficiency, i.e., the amount of energy spent to execute one single operation. Computational efficiency in classical computers levels off below \(10\) GMAC/s/W [6].

An alternative computing architecture with a more distributed interconnectivity and memory would allow for greater energy efficiency and computational speed. An inspiring example would be an architecture such as the brain. The brain is able to perform about \(10^{18}\) MAC/s using only \(20\,\mathrm{W}\) of power [6], and operates approximately \(10^{11}\) neurons with an average number of inputs for each of about \(10^{4}\) synapses. This leads to an estimated total of \(10^{15}\) synaptic connections, all conveying signals up to \(1\,\mathrm{kHz}\) bandwidth. The brain's computational efficiency (being less than \(1\,\mathrm{aJ}\) per MAC) is then about \(8\) orders of magnitude higher than the one of current supercomputers, which operate instead at \(100\,\mathrm{pJ}\) per MAC [6].

Connecting software to hardware through computing architecture tailored to ML tasks is the endeavor of research within the field of neuromorphic computing. The electronics community is now busy developing non-von Neumann computing architectures to enable information processing with an energy efficiency down to a few pJ per operation. Aiming to replicate fundamentals of biological neural circuits in dedicated hardware, important advances have been made in neuromorphic accelerators [7]. These advances are based on the spiking architectural models, which are still not fully understood. Deep Learning (DL)-focused approaches, on the other hand, aim to construct hardware that efficiently realizes DL architectures, while eliminating as much of the complexity of biological neural networks as possible. Among the most powerful DL hardware we can name the GPU-based DL accelerators hardware [8; 9; 10; 11; 12], as well as emerging analogue electronic Artificial Intelligence chipsets that tend to collocate processing and memory to minimize the memory-processor communication energy costs (e.g. the analogue
crossbar approaches [13]). The Mythic's architecture, for example, can yield high accuracy in inference applications within a remarkable energy efficiency of just half a pJ per MAC. Even if the implementation of neuromorphic approaches is visibly bringing outstanding record energy efficiencies and computation speeds, neuromorphic electronics is already struggling to offer the desired data throughput at the neuron level. Neuromorphic processing for high-bandwidth applications requires GHz operation per neuron, which calls for a fundamentally different technology approach.

### Optical Neural Networks

A major concern with neuromorphic electronics is that the distributed hardware needed for parallel interconnections is impractical to realize with classical metal wiring: a trade-off applies between interconnectivity and bandwidth, limiting these engine's utilization to applications in the kHz and sub-GHz regime. When sending information not through electrical signals but via optical signals, the optical interconnections do not undergo interference and the optical bandwidth is virtually unlimited. This can for example be achieved when exploiting the color and/or the space and/or the polarization and/or the time domain, thus allowing for applications in the GHz regime. It has been theorized that photonic neuromorphic processors could operate ten thousand times faster while using less energy per computation [14; 15; 16; 17]. Photonics therefore seems to be a promising platform for advances in neuromorphic computing.

Implementations of weighted addition for Optical Neural Networks (ONNs) include Mach-Zehnder Interferometer-based Optical Interference Units [18], time-multiplexed and, coherent detection [19], free space systems using spatial light modulators [20] and Micro-Ring-Resonator-based weighting bank on silicone [21]. Furthermore, Indium phosphide-integrated optical cross-connect using Semiconductor Optical Amplifiers as single stage weight elements, as well as Semiconductor Optical Amplifier-based wavelength converters [22; 23; 24] have been demonstrated for allowing All-Optical (AO) Neural Networks (NNs). A comprehensive review of all the approaches used in integrated photonics can be found in [25].

Next to these promises, aspects like implementation of nonlinearities, access and storage of weights in on-chip memory, and noise sources in analog photonic implementations, all pose challenges in devising scalable photonic neuromorphic processors and accelerators. These challenges also occur when they are embedded within end-to-end systems. Fortunately, arbitrary scalability of these networks has been demonstrated, with a certain noise and accuracy. However, it would be useful to envision new architectures to reduce noise even more.

### Noise in ONNs

The types of noise in ONNs include thermal crosstalk [26], cumulative noise in optical communication links [27; 28] and noise deriving from applying an activation function [29].

In all these studies, the noise is considered to be approximated well by Additive White Gaussian Noise (AWGN).

For example, taking the studies [26; 28; 27; 29; 30] as starting point, the authors of [31] model an ONN as a communication channel with AWGN. We follow this assumption and will model an ONN as having been built up from interconnected nodes with noise in between them. This generic approach does not restrict us to any specific device that may be used in practice.

The model also applies to the two alternative designs of an AO implementation of a NN (see for example [32]) and the case of an optical/electrical/optical (O/E/O) NN [22]. In an AO NN, the activation function is applied by manipulating an incoming electromagnetic wave. Modulation (and the AWGN it causes) only occurs prior to entering an AO NN (or equivalently, in the first layer). For the remainder of the network the signal remains in the optical domain. Here, when applying the optical activation function a new source of noise is introduced as AWGN at the end of each layer. Using the O/E/O network architecture, the weighted addition is performed in the optical realm, but the light is captured soon after each layer, where it is converted into an electrical and digital signal and the activation function is applied via software on a computer. The operation on the computer can be assumed to be noiseless. However, since the result again needs to be modulated (to be able to act as input to the next layer), modulation noise is added. We can further abstract from the specifics of the AO and O/E/O design and see that in either implementation noise occurs at the same locations within the mathematical modeling, namely AWGN for weighted addition and afterwards AWGN from an optical activation function or from modulation, respectively. This means that we do not need to distinguish between the two design choices in our modeling; we only need to choose the corresponding AWGN term after activation.

The operation of a layer of a feed-forward NN can be modeled by multiplying a matrix \(W\) with an input vector \(x\) (a bias term \(b\) can be absorbed into the matrix-vector product and will therefore suppressed in notation here) and then applying an activation function \(f:\mathbb{R}\rightarrow\mathbb{R}\) element-wise to the result. Symbolically,

\[x\mapsto f(Wx).\]

Now, concretely, the noise model that we study is described by

\[x\mapsto f(Wx+\mathrm{Normal}(0,\Sigma_{\mathrm{w}}))+\mathrm{Normal}(0,\Sigma _{\mathrm{a}}), \tag{1}\]

for each hidden layer of the ONN. Here \(\mathrm{Normal}(0,\Sigma)\) denotes the multivariate normal distribution with mean vector \(0\) and covariance matrix \(\Sigma\). More specifically, \(\Sigma_{\mathrm{w}}\), \(\Sigma_{\mathrm{a}}\) and \(\Sigma_{\mathrm{m}}\) are the covariance matrices associated with weighted addition, application of the activation function, and modulation, respectively. Figure 1 gives a schematic representation of the noise model under study. As we have seen above, in the O/E/O case we have \(\Sigma_{\mathrm{a}}=\Sigma_{\mathrm{m}}\), otherwise \(\Sigma_{\mathrm{a}}\) is due to the specific structure of the photonic activation function. The first layer, regardless of an AO or O/E/O network, sees a modulated input \(x\), i.e., \(x+\mathrm{Normal}(0,\Sigma_{\mathrm{m}})\), and afterwards the same
steps of weighing and applying an activation function, that is (1). Arguably the hidden layers and their noise structure are the most important parts, especially in deep NNs. Therefore, the main equation governing the behavior of the noise propagation in an ONN will remain (1).

### Noise-resistant designs for ONNs

The main contribution of this paper lies in analyzing two noise reduction mechanisms for feed-forward ONNs. The mechanisms are derived from the insight that noise can be mitigated through averaging because of the law of large numbers, and they are aimed at using the enormous bandwidth that photonics offer. The first design (Design A) and its analysis are inspired by recent advancements for NNs with random edges in [33]; the second design (Design B) is new and simpler to implement, but comes without a theoretical guarantee of correctness for nonlinear ONNs, specifically.

Both designs--illustrated in Figure 2--are built from a given NN for which an optical implementation is desired. Each design proposes a larger ONN by taking parts of the original NN, and duplicating and arranging them in a certain way. If noise is absent, then this larger ONN produces the same output as the original NN; and, if noise is present, then this ONN produces an output closer to the desired NN than the direct implementation of the NN as an ONN without modifications would give.

The first mechanism to construct a larger ONN suppressing inherent noise of analog systems starts with a certain number of copies \(N\) of the input data. The copies are all processed independently by (in parallel arranged copies of) the layers. Each copy of a layer takes in multiple input copies to produce the result of weighted addition, to which the activation mechanism is applied. The copies that are transmitted to each layer (or set of parallel arrayed layers) are independent of each other. The independent outputs function as inputs to the upcoming (copies of the) layers, and so on and so forth.

The idea of the second design is to use multiple copies of the input, on which the weighted addition is performed. The noisy products of the weighted addition are averaged to a single number/light beam. This average is then copied and the multiple copies are fed through the activation function, creating multiple noisy activations to be used as the next layer's input, and so on.

### Summary of results

Using Design A, we are able to establish that ONNs posses the same theoretical properties as NNs. Specifically, we can prove that _any_ NN can be approximated arbitrarily well by an ONN built using Design A (Theorem 1). Similar considerations for NNs with random edges can be found in [33], but the noise model and proof method are different. Here, we first bound the deviation of an ONN and a noiseless NN. To this bound Hoeffding's inequality is then applied.

Establishing this theoretical guarantee, however, is done by increasing the number of components exponentially as the depth of the network increases. The current proof shows that for an ONN with Design A meant to approximate a NN with \(L\) layers arbitrarily well (and thus reduce the noise to negligible levels), a sufficient number of components is \(\omega(K^{L(L+1)}L^{L})\) for some constant \(K>0\). This is however not to say that such a large number is necessary: it is merely sufficient.

From a practical viewpoint, however, having to use as few components as possible would be more attractive. We therefore also investigate Design B, in which the number of components increases only linearly with the depth of the network. Because Design A already allows us to establish the approximation property of ONNs, we limit our analysis of Design B to _linear_ NNs for simplicity. We specifically establish in Theorem 2 for any linear NN the exact output distribution of an ONN built using Design B. Similar to the guarantee for Design A in Theorem 1, but more restrictively, this implies that any linear NN can be approximated arbitrarily well by some ONN built using Design B. Strictly speaking, Design B now has no guarantee of correctness for nonlinear NNs, but this should practically not withhold us (especially when activations, for instance, are close to linear).

We conduct numerical experiments with Designs A and B by constructing to LeNet ONNs. The numerical results indicate that in practice, adding some components for noise negation is already sufficient to increase the accuracy of an

Figure 1: Schematic depiction of the noise model of ONNs that we study. First, data \(x\) is modulated onto light. This step adds an AWGN term \(N_{\text{m}}\). This light enters the Photonic Layer, in which a weighted addition takes place, adding AWGN \(N_{\text{w}}\). The activation function is then applied, adding AWGN \(N_{\text{a}}\). The activation function may be applied by photo-detecting the signal of the weighted addition, turning it to a digital signal and applying the activation function on a computer. The result of that action would then be modulated again, to produce the optical output of the photonic neuron. The modulator is thus only required in the first layer, as each photonic neuron takes in light and outputs light.
ONN; an exponential number does not appear not to be necessary (see Figures 3 to 4).

Finally, we want to remark that the high bandwidth of photonic circuits can be exploited to implement the designs as efficiently as possible.

### Outline of the paper

We introduce the AWGN model formally in Section 2. This model is the basis for the analysis of the proposed noise reduction schemes that are next discussed in Sections 3 and 4. There, we specifically define Designs A and B, and each design is followed by a mathematical analysis. The main results are Theorems 1 and 2. Section 5 contain numerical simulations on LeNet ONNs to which we apply Designs A and B. Section 6 concludes; technical details are deferred to the Appendix.

## 2 Model

We consider general feed-forward NNs implemented on analog optical devices. Noise occurs due to various reasons in those optical devices. Reasons include quantum noise in modulation, chip imperfections, and crosstalk [26; 28; 27; 29; 30].

The noise profiles and levels of different devices differ, but we can, to good approximation, expect AWGN to occur at three separate instances [31]: when modulating, when weighting, and when applying an activation function. The thus proposed AWGN model is formalized next in Section 2.1.

### Feed-forward nonlinear ONNs

We assume that our aim is to implement a feed-forward nonlinear NN with domain \(\mathbb{R}^{d_{0}}\) and range \(\mathbb{R}^{d_{L}}\), that can be represented by a parameterized function \(\Psi^{\text{NN}}:\mathbb{R}^{d_{0}}\times\mathbb{R}^{n}\to\mathbb{R}^{d_{L}}\) as follows. For \(\ell=1,\ldots,L\in\mathbb{N}_{+}\), \(\Psi^{\text{NN}}\) must be the composition of the functions

\[\Psi^{\text{NN}}_{\ell}:\mathbb{R}^{d_{\ell-1}}\to\mathbb{R}^{d_{\ell}},\quad x \mapsto\sigma^{(\ell)}\big{(}W^{(\ell)}x+b^{(\ell)}\big{)}.\]

Here, \(W^{(\ell)}\in\mathbb{R}^{d_{\ell}\times d_{\ell-1}}\) denotes the weight matrix in the \(\ell\)-th layer, \(b^{(\ell)}\in\mathbb{R}^{d_{\ell}\times 1}\) the bias vector in the \(\ell\)-th layer, and \(\sigma^{(\ell)}:\mathbb{R}^{d_{\ell}\times 1}\to\mathbb{R}^{d_{\ell}\times 1}\) the activation function in the \(\ell\)-th layer. Specifically, the NN satisfies

\[\Psi^{\text{NN}}(\cdot\,,w)=\Psi^{\text{NN}}_{L}(\cdot\,,w^{(L)})\circ\cdots \circ\Psi^{\text{NN}}_{1}(\cdot\,,w^{(1)}), \tag{2}\]

where \(w^{(\ell)}=(W^{(\ell)},b^{(\ell)})\) represents the parameters in the \(\ell\)-th layer. Note that we do not necessarily assume that the activation function is applied component-wise (it could be any high-dimensional function). Such cases are simply contained within the model.

Suppose now that the NN in (2) is implemented as an ONN, but _without_ amending its design. AWGN will then disrupt the output of each layer. Specifically, for depths \(L\in\mathbb{N}_{+}\), the ONN will be representable by a function \(\Psi^{\text{ONN}}\) that is the composition of the noisy functions

\[\Psi^{\text{ONN}}_{\ell}:\mathbb{R}^{d_{\ell-1}}\to\mathbb{R}^{d_{\ell}}, \quad x\mapsto\sigma^{(\ell)}\big{(}W^{(\ell)}x+b^{(\ell)}+N^{(\ell)}_{\text{ w}}\big{)}+N^{(\ell)}_{\text{a}} \tag{3}\]

Figure 2: (a) Base \(4-3-2\) network, light circles indicate activations, boxes indicate post-activations. (b) Example for Design A with \(2\) layers as input copies to each subsequent layer. The light circles indicate the linear operations/matrix-vector products. The results of the linear operation is averaged (single solid-blue circle) and fed through the activation function, producing the multiple version of the layers output (boxes). (c) Example of Design B.
for \(\ell=1,\ldots,L\in\mathbb{N}_{+}\). Here,

\[N_{\mathsf{w}}^{(\ell)}\stackrel{{\mathsf{(d)}}}{{=}}\mathrm{Normal}( 0,\Sigma_{\mathsf{w}}^{(\ell)})\text{ and }N_{\mathsf{a}}^{(\ell)}\stackrel{{ \mathsf{(d)}}}{{=}}\mathrm{Normal}(0,\Sigma_{\mathsf{act}}^{(\ell)})\]

denote multivariate normal distributions that describe the AWGN within the ONN. In other words, the ONN will satisfy

\[\Psi^{\mathrm{ONN}}(\cdot\,,w)=\Psi^{\mathrm{ONN}}_{L}(\cdot\,,w^{(L)})\circ \cdots\circ\Psi^{\mathrm{ONN}}_{1}(\cdot\,,w^{(1)}) \tag{4}\]

instead of (2). Observe that (4) is a random NN; its outcome is uncertain, but hopefully close to that of (2).

### Feed-forward linear ONNs

Let us briefly examine the special case of a feed-forward _linear_ ONN in more detail. That is, we now assume additionally that for \(\ell=1,\ldots,L\), there exist \(e^{(\ell)}\in\mathbb{R}^{d_{\ell}}\) such that \(\sigma^{(\ell)}(y)=D^{(\ell)}y\) where \(D^{(\ell)}=\mathrm{diag}(e^{(\ell)})\). In other words, each activation function \(\sigma^{(\ell)}\) does element-wise multiplications by constants.

If each activation function is linear, then the output distribution of each layer will remain multivariate normal distributed due to the so-called linear transformation theorem (34, Theorem 1.2.6). The mean and covariance matrix of the underlying multivariate normal distribution will however be transformed in each layer.

Let us illustrate how the covariance matrix transforms by discussing the first layer in detail. Each layer in (3) can be interpreted as a random function that takes the noisy vector \(\mathbf{A}^{(\ell-1)}=(\mathbf{A}_{1}^{(\ell-1)},\ldots,\mathbf{A}_{d_{\ell-1}} ^{(\ell-1)})\) say as input, and produces the even noisier vector \(\mathbf{A}^{(\ell)}=(\mathbf{A}_{1}^{(\ell)},\ldots,\mathbf{A}_{d_{\ell}}^{( \ell)})\) say as output. Specifically, the noisy input to the first layer is modeled by

\[\mathbf{A}^{(0)}\mid x\stackrel{{\mathsf{(d)}}}{{=}}x+\mathscr{N} \left(0,\Sigma_{\mathrm{m}}\right) \tag{5}\]

because of the modulation error within the first layer. Here \(\cdot\mid\cdot\) indicates a conditional random variable. This input next experiences weighted addition and more noise is introduced: the noisy preactivation of the first layer satisfies

\[\mathbf{U}^{(1)}\mid\mathbf{A}^{(0)}\stackrel{{ \mathsf{(d)}}}{{:=}}W^{(1)}\mathbf{A}^{(0)}+b^{(1)}+\mathscr{N}(0,\Sigma_{ \mathsf{w}}^{(1)}). \tag{6}\]

Combining (5) and (6) with the linear transformation theorem for the multivariate normal distribution as well as the fact that sums of independent multivariate normal random variables are again multivariate normally distributed (34, Theorem 1.2.14), we find that

\[\mathbf{U}^{(1)}\mid x \stackrel{{\mathsf{(d)}}}{{=}}W^{(1)}x+b^{(1)}+W^{( 1)}\mathscr{N}\big{(}0,\Sigma_{\mathsf{m}}\big{)}+\mathscr{N}\big{(}0,\Sigma_ {\mathsf{w}}^{(1)}\big{)}\] \[\stackrel{{\mathsf{(d)}}}{{=}}W^{(1)}x+b^{(1)}+ \mathscr{N}\big{(}0,W^{(1)}\Sigma_{\mathrm{m}}(W^{(1)})^{\intercal}+\Sigma_{ \mathsf{w}}^{(1)}\big{)}.\]

After applying the linear activation function, we obtain

\[\mathbf{A}^{(1)}\mid x\stackrel{{\mathsf{(d)}}}{{=}} \sigma_{1}(\mathbf{U}^{(1)})+\mathscr{N}(0,\Sigma_{\mathsf{a}}^{(1)})\mid x\] \[\stackrel{{\mathsf{(d)}}}{{=}}D^{(1)}(W^{(1)}x+b^{(1)})\] \[\quad+\mathscr{N}\Big{(}0,\Sigma_{\mathsf{a}}^{(1)}+D^{(1)}\big{(} W^{(1)}\Sigma_{\mathrm{m}}(W^{(1)})^{\intercal}+\Sigma_{\mathsf{w}}^{(1)} \big{)}(D^{(1)})^{\intercal}\Big{)}\] \[\quad=\mathscr{N}(\Psi^{\mathrm{NN}}_{1}(x,w),\Sigma_{\mathrm{ONN }}^{(1)})\]

say. Observe that the unperturbed network's output remains intact, and is accompanied by a centered normal distribution with an increasingly involved covariance matrix:

\[\Sigma_{\mathrm{ONN}}^{(1)}\] \[=D^{(1)}\big{(}W^{(1)}\Sigma_{\mathrm{m}}(W^{(1)})^{\intercal}+ \Sigma_{\mathsf{w}}^{(1)}\big{)}(D^{(1)})^{\intercal}+\Sigma_{\mathsf{a}}^{(1)}\] \[=D^{(1)}W^{(1)}\Sigma_{\mathrm{m}}(W^{(1)})^{\intercal}(D^{(1)})^ {\intercal}+D^{(1)}\Sigma_{\mathsf{w}}^{(1)}(D^{(1)})^{\intercal}\] \[\quad+\Sigma_{\mathsf{a}}^{(1)}. \tag{7}\]

Observe furthermore that the covariance matrix in (7) is independent of the bias \(b^{(1)}\).

The calculations in eqs. (5) to (7) can readily be extended into a recursive proof that establishes the covariance matrix of the entire linear ONN. Specifically, for \(\ell=1,\ldots,L\), define the maps

\[T^{(\ell)}(\Sigma) :=D^{(\ell)}W^{(\ell)}\Sigma(W^{(\ell)})^{\intercal}(D^{(\ell)})^{\intercal}\] \[\quad+D^{(\ell)}\Sigma_{\mathsf{w}}^{(\ell)}(D^{(\ell)})^{ \intercal}+\Sigma_{a}^{(\ell)}. \tag{8}\]

We then have the following:

**Proposition 1** (Distribution of linear ONNs): _Assume that there exist vectors \(e^{(\ell)}\in\mathbb{R}^{d_{\ell}}\) such that \(\sigma^{(\ell)}(y)=\text{diag}(e^{(\ell)})y\). The feed-forward linear ONN in (4) then satisfies_

\[\Psi^{\mathrm{ONN}}(\cdot,w)\stackrel{{\mathsf{(d)}}}{{=}} \mathscr{N}\big{(}\Psi^{\mathrm{NN}}(\cdot,w),\Sigma_{\mathrm{ONN}}^{(L)} \big{)},\]

_where for \(\ell=L,L-1,\ldots,1\),_

\[\Sigma_{\mathrm{ONN}}^{(\ell)}=T^{(\ell)}(\Sigma_{\mathrm{ONN}}^{(\ell-1)}); \quad\text{and}\quad\Sigma_{\mathrm{ONN}}^{(0)}=\Sigma_{\mathrm{m}}.\]

In linear ONNs with symmetric noise (that is, the AWGN of each layer's noise sources has the same covariance matrix), Proposition 1's recursion simplifies. Introduce \(P^{(\ell)}:=\prod_{i=\ell+1}^{L}D^{(i)}W^{(i)}\) for notational convenience. The following is proved in Appendix A.1.1:

**Corollary 1** (Symmetric noise case): _Within the setting of Proposition 1, assume additionally that for all \(\ell\in\mathbb{N}_{+}\), \(\Sigma_{\mathsf{a}}^{(\ell)}=\Sigma_{\mathsf{a}}\) and \(\Sigma_{\mathsf{w}}^{(\ell)}=\Sigma_{\mathsf{w}}\). Then,_

\[\Sigma_{\mathrm{ONN}}^{(L)} =P^{(0)}\Sigma_{\mathrm{m}}(P^{(0)})^{\intercal}+\sum_{\ell=1}^{L}P ^{(\ell)}\Sigma_{\mathsf{a}}(P^{(\ell)})^{\intercal}\] \[\quad+\sum_{\ell=1}^{L}P^{(\ell)}D^{(\ell)}\Sigma_{\mathsf{w}}(D^ {(\ell)})^{\intercal}\big{(}P^{(\ell)}\big{)}^{\intercal}.\]

_If moreover for all \(\ell\in\mathbb{N}_{+}\), \(W^{(\ell)}=W\), \(D^{(\ell)}=D\), and \(\|D\|_{\mathrm{F}}\|W\|_{\mathrm{F}}<1\), then_

\[\lim_{L\to\infty}\Sigma_{\mathrm{ONN}}^{(L)}=\sum_{n=0}^{\infty}(DW)^{n}\left(D \Sigma_{\mathsf{w}}D^{\intercal}+\Sigma_{\mathsf{a}}\right)\left((DW)^{n} \right)^{\intercal}.\]

Proposition 1 and Corollary 1 describe the output distribution of linear ONNs completely.
### Discussion

One way to think of the AWGN model in Section 2.1 is to take a step back from the microscopic analysis of individual devices, and consider an ONN as a series of black box devices (recall also Figure 1). Each black box device performs their designated task and acts as communication channels with AWGN. This way of modeling in order to analyze the impact of noise can also be seen in [31]; and other papers modeling optical channels include [28; 27]. Further papers considering noise in optical systems with similar noise assumptions are [35; 36], where furthermore multiplicative noise is considered when an amplifier is present within the circuit [35]. Qualitatively the results for Design A also apply for multiplicative noise, the scaling however may differ.

Limitations of the modelWe note firstly that modeling the noise in ONNs as AWGN is warranted only in an operating regime with many photons, and is thus unlikely to be a good model for ONNs that operate in a regime with just a few photons.

Secondly, due to physical device features and operation conditions, weights, activations, and outputs can only be realized in ONNs if their values lie in certain ranges. Such constraints are no part of the model in Section 2. Fortunately, however, the implied range restrictions are usually not a problem in practice. For example, activation functions like sigmoid and \(\tanh\) map into \([0,1]\) and \([-1,1]\), respectively. Additional regularization rules like weight decay also move the entries of weight matrices in NNs towards smaller values. In case physical constraints were met one can increase the weight decay parameter to further penalize large weights during training, leading to smaller weights so that the ONN is again applicable.

## 3 Results--Design A

### Reducing the noise in feed-forward ONNs (Design A)

Recall that an example of Design A is presented in Figure 2(b). Algorithm 1 constructs this tree-like network, given the desired number of copies \(n_{0},\ldots,n_{L}\) per layer.

```
0:Input \(\mathbf{n}=(n_{\ell})_{\ell=0,\ldots,L}\)
0:\(\prod_{\ell=0}^{L}n_{i}\) copies of input \(x^{(0)}\), named \({}^{1}x^{(0)},\ldots,\big{(}\prod_{\ell=0}^{L}n_{i}\big{)}\,x^{(0)}\) for\(\ell=0,\ldots,L-1\)do for\(\alpha=1,\ldots,\prod_{\ell=1}^{L-1}n_{i}\)do \({}^{\alpha}\xi^{(\ell)}\gets W^{(\ell+1)}{}^{\alpha}x^{(\ell)}+b^{(\ell+1 )}+\mathrm{Normal}(0,\Sigma_{\mathbf{w}})\) endfor for\(\alpha=0,\ldots,\big{(}\prod_{i=\ell}^{L-1}n_{i}\big{)}-1\)do \({}^{\alpha}y^{(\ell)}\xrightarrow{\mathrm{averaging}}n_{i}^{-1}\left({}^{ \alpha-n_{\ell}+1}\xi^{(\ell)}+\cdots+{}^{\alpha-n_{\ell}+n_{\ell}}\,\xi^{( \ell)}\right)\) \({}^{\alpha}x^{(\ell+1)}\leftarrow\sigma^{(\ell)}({}^{\alpha}y^{(\ell)})+ \mathrm{Normal}(0,\Sigma_{\mathbf{a}})\) endfor endfor return\({}^{1}x^{(L)}\)
```

**Algorithm 1** Algorithm to construct a noise reducing network

Observe that in Design A, the number of copies utilized in each layer, the \(n_{\ell}\), are fixed. There is however only a single copy in the last layer. Its output is the unique output of the ONN. Each other layer receives multiple independent inputs. With each of the independent copies weighted addition is performed, and the results are averaged to produce the layer's single output. Having independent incoming copies is achieved by having multiple independent branches of the prior partial networks incoming into a given layer. This means that the single layer \(L\) receives \(n_{L-1}\) independent inputs of \(n_{L-1}\) independent layers \(L-1\). Each of the \(n_{L-1}\) copies of layer \(L-1\) receives \(n_{L-2}\) inputs from independent copies of layer \(L-2\). Generally, let \(n_{\ell-1}\) be the number of copies of layer \(\ell-1\) that act as inputs to layer \(\ell\).

Observe that all copies are created upfront. That means there are \(\prod_{\ell=0}^{L-1}n_{\ell}\) copies of the data. By Algorithm 1, \(\prod_{\ell=1}^{L-1}n_{\ell}\) copies of the first layer are arrayed in parallel to each other, and each of them processes \(n_{0}\) copies of the data. The outputs of the \(\prod_{\ell=1}^{L-1}n_{\ell}\) arrayed copies of the first layer are the input to the \(\prod_{\ell=2}^{L-1}n_{\ell}\) arrayed copies of the second layer, and so on.

Notice that noise stemming from applying the activation function is subject to a linear transformation in the next layer. The activation function noise can therefore be considered as weight-noise by inserting an identity layer with \(\sigma=\mathrm{id}\), \(W=I\) and \(b=0\).

We want to verify that a Design A ONN yields outputs that are with high probability close to the original noiseless NN. Let \(\tilde{\Psi}^{\mathrm{ONN}}(x,w)\) the Design A ONN and then let

\[\mathbb{P}\Big{[}\sup_{x\in\mathbb{R}^{d}}\big{\|}\Psi^{\mathrm{ NN}}(x,w)-\tilde{\Psi}^{\mathrm{ONN}}(x,w)\big{\|}_{2}<D_{L}\Big{]}>1-C_{L}, \tag{9}\]

be the desired property. The main result of this section is the following:

**Theorem 1**: _For any \(C_{L}\in(0,1)\), any \(D_{L}\in(0,\infty)\), and any nonlinear NN \(\Psi^{\mathrm{NN}}\), with Lipschitz-continuous activations functions with Lipschitz-constants \(a^{(i)}\) and weight-matrices \(W^{(i)}\), Algorithm 1 is able to construct an ONN \(\tilde{\Psi}^{\mathrm{ONN}}\) that satisfies (9)._

_Let the covariance matrices of the occurring AWGN be diagonal matrices and let each of the values of the covariance matrices be upper bounded by \(\sigma^{2}\geq 0\). For any set of \((\kappa_{i})_{i=1,\ldots,L}\), \((\delta_{i})_{i=1,\ldots,L}\) such that \(\prod(1-\kappa_{\ell})>1-C_{L}\) and \(\sum\delta_{\ell}\leq D_{L}\), a sufficient number of copies to construct an ONN \(\tilde{\Psi}^{\mathrm{ONN}}\) that satisfies (9) is given by_

\[n_{L} =1\] \[n_{\ell} \geq\frac{\sigma^{2}\Big{(}\prod_{i=\ell+1}^{L}a^{(i)}\prod_{i= \ell+2}^{L}\|W^{(i)}\|_{\mathrm{op}}\Big{)}^{2}}{\delta_{\ell+1}^{2}}\] \[\quad\times\left(\sqrt{2}\frac{\Gamma((d_{\ell+1}+1)/2)}{\Gamma(d_ {\ell+1}/2)}\right.\]
\[\leq a^{(L)}\Big{\|}\frac{1}{n_{L-1}}\sum_{i=1}^{n_{L-1}}\] \[\quad\times\Big{\|}\sum_{i=1}^{n_{L-1}}\Big{(}\tilde{x}^{i}-\sigma^{ (L-1)}\big{(}W^{(L-1)}\big{(}\dots\big{)}+b^{(L-1)}\big{)}\Big{)}\Big{\|}_{2}\] \[\quad+a^{(L)}\Big{\|}\frac{1}{n_{L-1}}\sum_{i=1}^{n_{L-1}}N^{(i)} \Big{\|}_{2}.\]
In the next iteration step the term

\[\Big{\|}\sum_{i=1}^{n_{L-1}}\Big{(}\tilde{x}^{i}-\sigma^{(L-1)}\big{(}W^{(L-1)} \big{(}\dots\big{)}+b^{(L-1)}\big{)}\Big{)}\Big{\|}_{2}\]

is further bounded by first using the triangle inequality and thereafter bounding in the same way as we did in the first layer:

\[\Big{\|}\sum_{i=1}^{n_{L-1}}\Big{(}\sigma^{(L-1)}\Big{(}\frac{1}{n _{L-2}}\sum_{j_{i}=1}^{n_{L-2}}\Big{(}W^{(L-1)}\tilde{x}^{j_{i}}+b^{(L-1)}+N^{ (j_{i})}\Big{)}\Big{)}\] \[-\sigma^{(L-1)}\big{(}W^{(L-1)}\big{(}\sigma^{(L-2)}\big{(}W^{(L- 2)}\big{(}\dots\big{)}+b^{(L-2)}\big{)}\big{)}\] \[\qquad\qquad+b^{(L-1)}\big{)}\Big{)}\Big{\|}_{2}\] \[\leq\frac{a^{(L-1)}\|W^{(L-1)}\|_{\text{op}}}{n_{L-2}}\] \[\quad\times\sum_{i=1}^{n_{L-1}}\Big{\|}\sum_{j_{i}=1}^{n_{L-2}} \Big{(}\tilde{x}^{j_{i}}-\sigma^{(L-2)}\big{(}W^{(L-2)}\big{(}\dots\big{)}+b^ {(L-2)}\big{)}\Big{)}\Big{\|}_{2}\] \[\quad+a^{(L-1)}\sum_{i=1}^{n_{L-1}}\Big{\|}\frac{1}{n_{L-2}}\sum_ {j_{i}=1}^{n_{L-2}}N^{(j_{i})}\Big{\|}_{2}.\]

Here,

\[\sum_{i=1}^{n_{L-1}}\Big{\|}\sum_{j_{i}=1}^{n_{L-2}}\Big{(}\tilde{x}^{j_{i}}- \sigma^{(L-2)}\big{(}W^{(L-2)}\big{(}\dots\big{)}+b^{(L-2)}\big{)}\Big{)} \Big{\|}_{2}\]

may again be bounded in the same fashion. This leads to the following recursive argument.

Let \(\mathscr{F}^{(\ell)}\) be the sum of the differences between-- loosely speaking--the ends of the remaining Design A "subtrees" and noiseless NNs "subtrees" at layer \(\ell\). More specifically, let

\[\mathscr{F}^{(L)} :=\Big{\|}\tilde{x}-\sigma^{(L)}\big{(}W^{(L)}\big{(}\dots\big{)} +b^{(L)}\big{)}\Big{\|};\] \[\mathscr{F}^{(\ell)} :=\sum_{i_{L}=1}^{n_{L-1}}\sum_{i_{L-L}=1}^{n_{L-2}}\dots\sum_{i _{\ell+2}\dots i_{L-1}=1}^{n_{\ell+1}}\Big{\|}\sum_{i_{\ell+1\ell+2\dots i_{ L-1}=1}}^{n_{\ell}}\] \[\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad \qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad \qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad \qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad \qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad \qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad \qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad \qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad \qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad \qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad \qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad \qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad \qquad\
\[\mathbb{P}\left[\sum_{\ell}\mathscr{S}_{\ell}<D_{L}\right]\geq \mathbb{P}\left[\bigcap_{\ell}\{\mathscr{S}_{\ell}<\delta_{\ell}\}\right]\] \[=\prod_{\ell}\mathbb{P}\Big{[}\mathscr{S}_{\ell}<\delta_{\ell} \Big{]}>\prod_{\ell}(1-\kappa_{\ell})>1-C_{L}.\]

Here, in the first inequality the dependence on \(x\) disappears due to (14).

#### 3.3.2 Bound for deviations

We next consider the \(\mathscr{S}_{\ell}\) for which we want to guarantee that

\[\mathbb{P}\big{[}\mathscr{S}_{\ell}<\delta_{\ell}\big{]}>1- \kappa_{\ell}.\]

Let \(m_{\ell}=\prod_{k=\ell}^{L}n_{k}\). By assumption the \(N_{k}^{(j_{i})}\) are independent and identically \(\mathrm{Normal}(0,\sigma_{k}^{2})\) distributed, where \(\sigma_{k}^{2}\leq\sigma^{2}\), for some common \(\sigma^{2}\). We are lower bounding the number of copies required, therefore using AWGN with higher variance only increases the lower bound, as the calculations below show. We calculate the bound exemplary for \(N^{(j_{i})}\) distributed according to \(\mathrm{Normal}(0,\sigma^{2})\), re-substituting \(\sigma_{k}^{2}\) below in (18) (which is the bound given in Theorem 1) thus covers the case of \(N_{k}^{(j_{i})}\stackrel{{\text{\sf(d)}}}{{=}}\mathrm{Normal}(0, \sigma_{k}^{2})\).

Each component of the vector

\[\sum_{j_{i}=1}^{n_{\ell-1}}N^{(j_{i})}=\Big{(}\sum_{j_{i}=1}^{n_{ \ell-1}}N_{1}^{(j_{i})},\ldots,\sum_{j_{i}=1}^{n_{\ell-1}}N_{d_{\ell}}^{(j_{i} )}\Big{)}^{\intercal}\]

is assumed to be \(\mathrm{Normal}(0,n_{\ell-1}\sigma^{2})=\sqrt{n_{\ell-1}}\sigma\mathrm{Normal }(0,1)\) distributed. It then holds that

\[\sum_{i=1}^{m_{\ell}}\Bigl{\|}\sum_{j_{i}=1}^{n_{\ell-1}}N^{(j_{i })}\Bigr{\|}_{2}\stackrel{{\text{\sf(d)}}}{{=}}\sum_{i=1}^{m_{ \ell}}\sqrt{n_{\ell-1}}\sigma\|\mathrm{Normal}(0,I_{d})\|_{2}.\]

This is a sum of independent chi-distributed random variables, which means they are sub-gaussian (see below that we can calculate the sub-gaussian norm and it is indeed finite). Thus Hoeffding's inequality applies, according to which, for \(X_{1},\ldots,X_{n}\) independent, mean zero, sub-gaussian random variables, for every \(t\geq 0\)

\[\mathbb{P}\Big{[}\Bigl{|}\sum_{i=1}^{N}X_{i}\Bigr{|}<t\Bigr{]}>1- 2\exp\Bigl{(}-\frac{ct^{2}}{\sum_{i=1}^{N}\|X_{i}\|_{\psi_{2}}^{2}}\Bigr{)} \tag{16}\]

holds; see e.g. (37, Theorem 2.6.2). Here \(c>0\) is an absolute constant (see (37, Theorem 2.6.2)) and

\[\|X\|_{\psi_{2}}:=\inf\bigl{\{}t>0:\mathbb{E}[\exp(X^{2}/t^{2})] \leq 2\bigr{\}}.\]

To apply Hoeffding's inequality in our setting, we need to center the occurring random variables. For \(N^{(i)}\sim\mathrm{Normal}(0,I_{d})\), the term \(\|N^{(i)}\|_{2}\) is chi distributed with mean

\[\mu_{d}=\sqrt{2}\frac{\Gamma((d+1)/2)}{\Gamma(d/2)}, \tag{17}\]

where \(\Gamma\) is the gamma function, see e.g. (38, p.238).

Consider

\[\mathbb{P}\left[\frac{\prod_{i=\ell}^{L}a^{(i)}\prod_{i=\ell+1}^{ L}\|W^{(i)}\|_{\mathrm{op}}}{m_{\ell}\sqrt{n_{\ell-1}}}\sigma\sum_{i=1}^{m_{ \ell}}\|N^{(i)}\|_{2}<\delta_{\ell}\right]\] \[=\mathbb{P}\left[\frac{\prod_{i=\ell}^{L}a^{(i)}\prod_{i=\ell+1}^ {L}\|W^{(i)}\|_{\mathrm{op}}}{m_{\ell}\sqrt{n_{\ell-1}}}\sigma\sum_{i=1}^{m_{ \ell}}\Bigl{(}\|N^{(i)}\|_{2}-\mu_{d_{\ell}}\Bigr{)}\right.\] \[\left.\hskip 56.905512pt<\delta_{\ell}-\frac{\prod_{i=\ell}^{L}a^{( i)}\prod_{i=\ell+1}^{L}\|W^{(i)}\|_{\mathrm{op}}}{m_{\ell}\sqrt{n_{\ell-1}}} \sigma m_{\ell}\mu_{d_{\ell}}\right]\]

which equals

\[\mathbb{P}\Bigg{[}\sum_{i=1}^{m_{\ell}} \Big{(}\|N^{(i)}\|_{2}-\mu\Big{)}\] \[<\frac{m_{\ell}\sqrt{n_{\ell-1}}\delta_{\ell}}{\sigma\prod_{i= \ell}^{L}a^{(i)}\prod_{i=\ell+1}^{L}\|W^{(i)}\|_{\mathrm{op}}}-m_{\ell}\mu \Bigg{]}\]

and is lower bounded (compare to (16)) by

\[1-2\exp\left(\frac{-c\Big{(}\frac{m_{\ell}\sqrt{n_{\ell-1}}\delta_{\ell}}{ \sigma\prod_{i=\ell}^{L}a^{(i)}\prod_{i=\ell+1}^{L}\|W^{(i)}\|_{\mathrm{op}}} -m_{\ell}\mu_{d_{\ell}}\Big{)}^{2}}{\sum_{i=1}^{m_{\ell}}C^{2}\Bigl{\|}N^{(i)} \|_{2}-\mu\Big{\|}_{\psi_{2}}^{2}}\right),\]

which in turn is lower bounded by

\[1-2\exp\left(\frac{-c\Big{(}\frac{m_{\ell}\sqrt{n_{\ell-1}}\delta_{\ell}}{ \sigma\prod_{i=\ell}^{L}a^{(i)}\prod_{i=\ell+1}^{L}\|W^{(i)}\|_{\mathrm{op}}} -m_{\ell}\mu_{d_{\ell}}\Big{)}^{2}}{\sum_{i=1}^{m_{\ell}}C^{2}\Bigl{\|}N^{(i)} \|_{2}\Bigr{\|}_{\psi_{2}}^{2}}\right),\]

where \(C>0\) is an absolute constant (see (37, Lemma 2.6.8)). For a chi distributed random variable \(\mathbf{X}\) it holds that

\[\mathbb{E}[\exp(\mathbf{X}^{2}/t^{2})]=M_{\mathbf{X}^{2}}(1/t^{2})\]

where \(M_{\mathbf{X}^{2}}(s)\) is the moment generating function of \(\mathbf{X}^{2}\)--a chi-squared distributed random variable. It is known (see e.g. (39, Appendix 13)) that

\[M_{\mathbf{X}^{2}}(s)=(1-2s)^{-d_{\ell}/2}\]

for \(s<\frac{1}{2}\). Accordingly for \(2<t^{2}\), the property in the definition of the sub-gaussian norm

\[\mathbb{E}[\exp(\mathbf{X}^{2}/t^{2})]=\Big{(}1-2\frac{1}{t^{2}} \Big{)}^{-d_{\ell}/2}\leq 2\]

is satisfied for all \(t\) for which

\[t\geq\max\left\{\sqrt{\frac{4\sqrt[4]{4}}{2\sqrt[4]{4}-2}},\sqrt{2} \right\}=\sqrt{\frac{4\sqrt[4]{4}}{2\sqrt[4]{4}-2}}\]

holds. The square of the sub-gaussian norm of the chi distributed random variables is thus

\[\Bigl{\|}\|N^{(i)}\|_{2}\Bigr{\|}_{\psi_{2}}^{2}=\frac{4\sqrt[4]{4}}{2\sqrt[4] {4}-2}.\]
Substituting the norm into the lower bound yields

\[1-2\exp\left(\frac{-c\Big{(}\frac{m_{\ell}\sqrt{n_{\ell-1}}\delta_{\ell}}{\sigma \prod_{i=a^{(i)}}^{L}\prod_{i=\ell+1}^{\ell-1}\|W^{(i)}\|_{\text{op}}}-m_{\ell} \mu_{d_{\ell}}\Big{)}^{2}}{C^{2}m_{\ell}\frac{4\sqrt{4}}{2}}\right).\]

In order to achieve (15), a sufficient criterion is

\[\frac{\kappa_{\ell}}{2}\geq\exp\left(\frac{-cm_{\ell}\Big{(}\frac{\sqrt{n_{ \ell-1}}\delta_{\ell}}{\sigma\prod_{i=a^{(i)}}^{L}\prod_{i=\ell+1}^{\ell-1}\|W^ {(i)}\|_{\text{op}}}-\mu_{d_{\ell}}\Big{)}^{2}}{C^{2}\frac{4\sqrt{4}}{2}}\right).\]

Solving for \(n_{\ell-1}\) leads to

\[n_{\ell-1}\geq \frac{\sigma^{2}\Big{(}\prod_{i=\ell}^{L}a^{(i)}\prod_{i=\ell+1}^{ L}\|W^{(i)}\|_{\text{op}}\Big{)}^{2}}{\delta_{\ell}^{2}}\] \[\times\left(\sqrt{C^{2}\frac{4\sqrt{4}}{2\sqrt[4]{4}-2}(-\ln( \kappa_{\ell}/2))\frac{1}{cm_{\ell}}}+\mu_{d_{\ell}}\right)^{2}. \tag{18}\]

If we substitute the expression in (17) for \(\mu_{d_{\ell}}\), (18) becomes the bound as seen in Theorem 1.

### Conclusion

Within the context of the model described in Section 2, we have established that any feed-forward NN can be approximated arbitrarily well by ONNs constructed using Design A. This is Theorem 1 in essence.

This result has two consequences when it comes to the physical implementation of ONNs. On the one hand, it is guaranteed that the theoretical expressiveness of NNs can be retained in practice. On the other hand, Design A allows one to improve the accuracy of a noisy ONN to a desired level, and in fact bring the accuracy arbitrarily close to that of any state-of-the-art feed-forward noiseless NNs. Let us finally remark that the high bandwidth of photonic circuits may be of use when implementing Design A.

## 4 Results--Design B

### Reducing noise in feed-forward linear ONNs (Design B)

Recall that an example of Design B is presented in Figure 2(c). Algorithm 2 constructs this network, given a desired number of copies \(m\) in each layer.

```
0:\(m\) copies of input \({}^{1}x^{(0)},\dots,{}^{m}x^{(0)}\) for\(\ell=1,\dots,L\)do for\(\alpha=1,\dots,m\)do \({}^{\alpha}\xi^{(\ell)}\gets W^{(\ell)}\,{}^{\alpha}x^{(\ell-1)}+b^{( \ell)}+\text{Normal}(0,\Sigma_{\mathsf{w}})\) endfor \(y^{(\ell)}\xrightarrow[]{\text{combining}}\,{}^{1}\xi^{(\ell)}+\dots+{}^{m}\, \xi^{(\ell)}\) \(\binom{{}^{1}y^{(\ell)}\dots}{}^{m}y^{(\ell)}\xrightarrow[]{\text{splitting}}\,m^{-1}y^{(\ell)}\) for\(\alpha=1,\dots,m\)do \({}^{\alpha}x^{(\ell)}\leftarrow\sigma^{(\ell)}({}^{\alpha}y^{(\ell)})+\text{Normal}(0, \Sigma_{\mathsf{act}})\) endfor endfor return\(m^{-1}\sum_{\alpha=1}^{m}\,{}^{\alpha}x^{(L)}\)
```

**Algorithm 2** Algorithm to construct a noise reducing network

Calculating the output of a NN by using Design B first requires to fix a number \(m\). The input data \(x^{(0)}\) is then modulated \(m\) times, creating \(m\) noisy realizations of the input \(({}^{\alpha}x^{(0)})_{\alpha=1,\dots,m}\). The weighted addition step and the activation function of each layer are singled out and copied \(m\) times. Both the copies of the weighted addition step and of the activation function of each layer are arrayed parallel to each other and performed on the \(m\) inputs, resulting in \(m\) outputs. The \(m\) parallel outputs of the weighted addition are merged to a single output, and afterwards split into \(m\) pieces. The \(m\) pieces are each send to one of the \(m\) activation function mechanisms for processing. The resulting \(m\) activation values are the output of the layer. If it is the last layer, the \(m\) activation values are merged to produce the final output. These steps are formally described in Algorithm 2. A schematic representation of Design B can be seen in Figure 2(c).

```
0:\(m\) copies of input \({}^{1}x^{(0)},\dots,{}^{m}x^{(0)}\) for\(\ell=1,\dots,L\)do for\(\alpha=1,\dots,m\)do \({}^{\alpha}\xi^{(\ell)}\gets W^{(\ell)}\,{}^{\alpha}x^{(\ell-1)}+b^{(\ell)}+ \text{Normal}(0,\Sigma_{\mathsf{w}})\) endfor \(y^{(\ell)}\xrightarrow[]{\text{combining}}\,{}^{1}\xi^{(\ell)}+\dots+{}^{m}\, \xi^{(\ell)}\) \(\binom{{}^{1}y^{(\ell)}\dots}{}^{m}y^{(\ell)}\xrightarrow[]{\text{splitting}}\,m^{-1}y^{(\ell)}\) for\(\alpha=1,\dots,m\)do \({}^{\alpha}x^{(\ell)}\leftarrow\sigma^{(\ell)}({}^{\alpha}y^{(\ell)})+\text{Normal}(0, \Sigma_{\mathsf{act}})\) endfor endfor return\(m^{-1}\sum_{\alpha=1}^{m}\,{}^{\alpha}x^{(L)}\)
```

**Algorithm 3** Algorithm to construct a noise reducing network

### Analysis of Design B

We now consider the physical and mathematical consequences of Design B.

Observe that in Design B, the \(m\) weighted additions of the \(\ell\)-th layer's input \(x^{(l-1)}\) result in realizations \(({}^{\alpha}\xi^{(\ell)})_{\alpha=1,\dots,m}\) of \(W^{(\ell)}x^{(\ell-1)}+b^{(\ell)}+\text{Normal}(0,\Sigma_{\mathsf{w}})\). These realizations are then combined resulting in

\[mW^{(\ell)}x^{(\ell-1)}+mb^{(\ell)}+\text{Normal}(0,m\Sigma_{\mathsf{w}})+\text{ Normal}(0,\Sigma_{\mathsf{sum}}).\]

Splitting the signal again into \(m\) parts, each signal carries information following the distribution

\[W^{(\ell)}x^{(\ell-1)}+b^{(\ell)} +\text{Normal}(0,m^{-1}\Sigma_{\mathsf{w}}+m^{-2}\Sigma_{\mathsf{ sum}})\] \[+\text{Normal}(0,\Sigma_{\mathsf{spl}}).\]

The mean of the normal distribution therefore is the original networks pre-activation obtained from this input (that is without perturbations). The covariance matrix of the normal distribution is \(m^{-1}\Sigma_{\mathsf{w}}+m^{-2}\Sigma_{\mathsf{sum}}+\Sigma_{\mathsf{spl}}\). Each of those signals is fed through the mechanism applying the activation function, yielding \(m\) noisy versions of the output, distributed according to

\[x^{(\ell)}\mid x^{(\ell-1)}\] \[\stackrel{{(\mathsf{d})}}{{=}}\sigma^{(\ell)}\big{(}W^ {(\ell)}x^{(\ell-1)}+b^{(\ell)}\] \[\qquad\qquad+\text{Normal}(0,m^{-1}\Sigma_{\mathsf{w}}+m^{-2} \Sigma_{\mathsf{sum}}+\Sigma_{\mathsf{spl}})\big{)}\] \[+\text{Normal}(0,\Sigma_{\mathsf{act}}).\]

The effect of Design B is thus that \(T^{(\ell)}(\Sigma)\) in (8) is replaced by

\[T^{(\ell)}_{m}(\Sigma):=\frac{1}{m}D^{(\ell)}W^{(\ell)}\Sigma\big{(}D^{(\ell)}W^ {(\ell)}\big{)}^{\mathsf{T}}+\frac{1}{m}D^{(\ell)}\Sigma_{\mathsf{w}}\big{(}D^ {(\ell)}\big{)}^{\mathsf{T}}\]
\[\leq\frac{1}{m^{2}}D^{(\ell)}\Sigma_{\mathrm{sum}}\big{(}D^{(\ell)} \big{)}^{\intercal}+D^{(\ell)}\Sigma_{\mathrm{spl}}\big{(}D^{(\ell)}\big{)}^{ \intercal}+\Sigma_{\mathrm{a}};\]

see Appendix A.2.2. Observe also that \(\Sigma_{a}^{(\ell)}\) can be written as \((1/m)m\Sigma_{a}^{(\ell)}\). Therefore, if we substitute the matrix \(\bar{\Sigma}_{a}^{(\ell)}=m\Sigma_{a}^{(\ell)}\) for \(\Sigma_{a}^{(\ell)}\) in \(T^{(\ell)}(\Sigma)\), we can write

\[T_{m}^{(\ell)}(\Sigma) =m^{-1}T^{(\ell)}(\Sigma)+m^{-2}D^{(\ell)}\Sigma_{\mathrm{sum}} \big{(}D^{(\ell)}\big{)}^{\intercal}\] \[\quad+D^{(\ell)}\Sigma_{\mathrm{spl}}\big{(}D^{(\ell)}\big{)}^{ \intercal}.\]

We have the following analogs to Proposition1 and Corollary1:

**Theorem 2** (Distribution of Design B): _Assume that there exist vectors \(a^{(\ell)}\in\mathbb{R}^{d_{\ell}}\) such that \(\sigma^{(\ell)}(y)=\text{diag}(a^{(\ell)})y\). The feed-forward linear ONN constructed using Design B with \(m\) copies then satisfies_

\[\Psi_{m}^{\mathrm{ONN}}(\cdot,w)\stackrel{{(\ref{eq:DNN})}}{{=}} \mathrm{Normal}\big{(}\Psi^{\mathrm{NN}}(\cdot,w),\Sigma_{\mathrm{ONN},m}^{(L) }\big{)},\]

_where for \(\ell=L,L-1,\dots,1\),_

\[\Sigma_{\mathrm{ONN},m}^{(\ell)}=T_{m}^{(\ell)}(\Sigma_{\mathrm{ONN},m}^{(\ell -1)});\quad\text{and}\quad\Sigma_{\mathrm{ONN},m}^{(0)}=\Sigma_{\mathrm{m}}.\]

Under the assumption of symmetric noise, a similar simplification of the recursion in Theorem2, similar to that in Proposition1, is possible. Assume \(\Sigma_{\mathrm{sum}}=\Sigma_{\mathrm{spl}}=0\). Introduce again \(P^{(\ell)}:=\prod_{i=\ell+1}^{L}D^{(i)}W^{(i)}\) for notational convenience. The following is proved in AppendixA.2:

**Corollary 2** (Symmetric noise case): _Assume that for all \(\ell\in\mathbb{N}_{+}\), \(\Sigma_{\mathrm{a}}^{(\ell)}=\Sigma_{\mathrm{a}}\) and \(\Sigma_{\mathrm{w}}^{(\ell)}=\Sigma_{\mathrm{w}}\). Then,_

\[\Sigma_{\mathrm{ONN},m}^{(L)} =\sum_{\ell=1}^{L}(m^{-1})^{L-\ell}P^{(\ell)}m\Sigma_{\mathrm{a} }(P^{(\ell)})^{\intercal}\] \[\quad+\sum_{\ell=1}^{L}(m^{-1})^{L-\ell}P^{(\ell)}D^{(\ell)} \Sigma_{\mathrm{w}}(D^{(\ell)})^{\intercal}\big{(}P^{(\ell)}\big{)}^{\intercal}\] \[\quad+(m^{-L})P^{(0)}\Sigma_{\mathrm{m}}(P^{(0)})^{\intercal}.\]

We will next consider the limit of the covariance matrix in a large, symmetric linear ONN with Design B, that we can grow infinitely deep. Algorithm2 is namely able to guarantee boundedness of the covariance matrix in such deep ONN if the parameter \(m\) is chosen appropriately:

**Corollary 3**: _Consider a linear ONN with Design B and parameter \(m\), that has \(L\) layers, and that satisfies the following symmetry properties: for all \(\ell\in\{1,\dots,L\}\), \(W^{(\ell)}=W\), \(D^{(\ell)}=D\), \(\Sigma_{\mathrm{a}}^{(\ell)}=\Sigma_{\mathrm{a}}\) and \(\Sigma_{\mathrm{w}}^{(\ell)}=\Sigma_{\mathrm{w}}\). Then, if \(\|D\|_{\mathrm{F}}\|W\|_{\mathrm{F}}<\sqrt{m}\), the limit \(\lim_{L\to\infty}\Sigma_{\mathrm{ONN}}^{(L)}\) exists._

_Moreover,_

\[\lim_{L\to\infty}\Sigma_{\mathrm{ONN},m}^{(L)}\] \[=\sum_{n=0}^{\infty}m^{-(n+1)}(DW)^{n}\left(D\Sigma_{\mathrm{w}}D ^{\intercal}+m\Sigma_{\mathrm{a}}\right)\left((DW)^{n}\right)^{\intercal}.\]

Notice that the bound on the number of copies needed for the covariance matrix of an ONN to converge to a limit is independent of e.g. the Frobenius norms of the covariance matrices that describe the noise distributions. This is because, here, we are not interested in bounding the covariance matrix to a specific level; instead, we are merely interested in the existence of a limit.

### Discussion & Conclusion

Compared to Theorem2's recursive description of the covariance matrix in any linear ONN with Design B, Corollary2 provides a series that describes the covariance matrix in any linear, symmetric ONN with Design B. While the result holds more restrictively, it is more insightful. For example, it allows us to consider the limit of the covariance matrix in an extremely deep ONNs (see Corollary3). Corollary3 suggests that in deep ONNs with Design B, one should choose \(m\approx\lceil(\|D\|_{\mathrm{F}}\|W\|_{\mathrm{F}})^{2}\rceil\) in order to control the noise and not be too inefficient with the number of copies.

These results essentially mean that in a physical implementation of an increasingly deep and linear ONN, the covariance matrix can be reduced (and thus remain bounded) by applying Design B with multiple copies. The quality of the ONN's output increases as the number of copies in Design B (or Design A for that matter) is increased. Finally, it is worth mentioning that Design B could potentially be implemented such that it leverages the enormous bandwidth of optics.

## 5 Simulations

We investigate the improvements of output quality achieved by Designs A and B on a benchmark example: the convolutional neural network LeNet [40]. As measure for quality we consider the Mean Squared Error (MSE)

\[\mathrm{MSE}=\frac{1}{n}\sum_{i=1}^{n}\Big{(}\Psi^{\mathrm{ONN}}(x^{(i)},w)- \Psi^{\mathrm{NN}}(x^{(i)},w)\Big{)}^{2}\]

and the prediction accuracy

\[\frac{\#\{\text{correctly classified images}\}}{\#\{\text{images}\}}.\]

### Empirical variance

We extracted plausible values for \(\Sigma_{\mathrm{w}}\) and \(\Sigma_{\mathrm{a}}\) from the ONN implementation [41] of a 2-layer NN for classification of the Modified National Institute of Standards and Technology (MNIST) database [42]. In [41], the authors trained the NN on a classical computer and implemented the trained weights afterwards in an ONN. They then tuned noise (with the same noise model as in Section2 of this paper) into the noiseless computer model, assuming that \(\Sigma_{\mathrm{w}}=\mathrm{diag}(\sigma_{\mathrm{w}}^{2})\) and \(\Sigma_{\mathrm{a}}=\mathrm{diag}(\sigma_{\mathrm{a}}^{2})\). They found \(\sigma_{\mathrm{w}}\in[0.08,0.1]\cdot d\) and \(\sigma_{\mathrm{a}}\in[0.1,0.15]\cdot d\) to reach the same accuracy levels as the ONN, where \(d\) denotes the diameter of the range.
### LeNet ONN: Performance when implemented via Design A and B

Convolutional NNs can be regarded as feedforward NNs by stacking the (2D or 3D) images into column vectors and arranging the filters to a weight matrix. Thus Design A and B are well-defined for Convolutional NNs. We apply the designs to LeNet5 [40], which is trained for classifying the handwritten digits in the MNIST dataset [42]. The layers are:

1. 2D convolutional layer with kernel size \(5\), stride \(1\) and \(2\)-padding. Output has \(6\) channels of \(28\)x\(28\) pixel representations, with the activation function being \(\tanh\);
2. average pooling layer, pooling \(2\)x\(2\) block, the output therefore is \(14\)x\(14\);
3. 2D convolutional layer with kernel size \(5\), stride \(1\) and no padding. The output has \(16\) channels of \(10\)x\(10\) pixel representations and the activation function is \(\tanh\);
4. average pooling layer, pooling \(2\)x\(2\) block, the output therefore is \(5\)x\(5\);
5. 2D convolutional layer with kernel size \(5\), stride \(1\) and no padding. The output has \(120\) channels of \(1\) pixel representations and the activation function used is \(\tanh\);
6. flattening layer, which turns the \(120\) one-dimensional channels into one \(120\)-dimensional vector;
7. dense layer with \(84\) neurons and \(\tanh\) activation function;
8. dense layer with \(10\) neurons and \(\operatorname{softmax}\) activation function.

Figures 3 and 4 show the MSE and the prediction accuracy of Design A and B for an increasing number of copies, respectively.

For simplicity we set all individual _copies_\(n_{i}\)_per layer_\(i\) in Design A to equal \(m\), that is \(n_{i}=m\) for all \(i\). The total number of copies that Design A starts with then is \(m^{L}\). Here \(L\) is equal to \(7\). In Design B the number of copies is \(m\)_per layer_ and the total number of copies is \(mL\). In the case of one copy the designs A and B are identical to the original network, while we focus on the effect once the designs deviate from the original network (\(m\geq 2\)).

The axis in Figures 3 and 4 denote the number of copies _per layer_. Here, we scale the copies per layer for Design A linearly, because the total amount of copies for Design A grows exponentially and we scale the copies per layer for Design B exponentially, because the total number of copies for Design B grows linearly. This way the comparison is on equal terms.

Figure 3 displays the MSE seen for LeNet, depending on the amount of copies for each design. In the trade-off between additional resources needed for the additional copies against the diminishing benefits of adding further copies, we see that, for both measures MSE (Figure 3) and relative accuracy (Figure 4), already \(2\) to \(5\) copies per layer yield good results. The relative accuracy in Figure 4 is scaled such that \(0\) corresponds to the accuracy of the original NN with noise profile (i.e., the ONN without modifications, we call this the original ONN) and \(1\) to the accuracy of the original NN without noise. The designs do not alter the fundamental operation of the original NN, therefore there should be no performance gain and the original NN's accuracy should be considered the highest achievable, thus constituting the upper bound in relative accuracy of \(1\). Likewise the lowest accuracy should be given by the original ONN, as there is no noise reduction involved.

### Effect of additional layers in LeNet

In order to investigate how the depth affects the noise at the output, while keeping the operation of the network the same to ensure the results are commensurable, we insert additional layers with identity matrix and identity activation function (we will call them identity layers) into a network.

Figure 4: Relative accuracy for Design A (top) and Design B (bottom) as function of copies on LeNet5 trained for MNIST classification. The pale area contains the \(56.5\)%-confidence intervals.

Figure 3: MSE(\(\cdot 10^{2}\)) for Design A (top) and Design B (bottom) as function of copies on LeNet5 trained for MNIST classification. The pale area contains the 95%-confidence intervals.
Specifically, we take networks with the LeNet architecture as in Section 5.2, using different activation functions, while fixing the output layer to be \(\operatorname{softmax}\). We then insert identity layers between layers 1 and 2, 3 and 4, 5 and 6, as well as between layers 6 and 7. For a fixed total of additional layers, the layers are inserted in the four spots between layers \(1\&2\), \(3\&4\), \(5\&6\), and \(7\&8\) according to the tuple

\[n\mapsto\left(\Big{\lfloor}\frac{n+3}{4}\Big{\rfloor},\Big{\lfloor}\frac{n+2} {4}\Big{\rfloor},\Big{\lfloor}\frac{n+1}{4}\Big{\rfloor},\Big{\lfloor}\frac{ n}{4}\Big{\rfloor}\right).\]

The insertion pattern is illustrated in Table 1:

Finally, we tune the variance terms of the covariance matrix in our noise model. The results are displayed in Figure 5.

In Figure 5, we observe that the \(\tanh\) and the \(\operatorname{ReLU}\) networks perform as expected. Additional noisy layers decrease the accuracy and thus the same level of performance can only be achieved if the variance is lower. This trend can also be seen in the linear network, but to a lesser extend.

### Simulations on effective values for Design B

According to Corollary 3, the covariance matrix of a linear \(\mathsf{ONN}\) constructed by Design B is bounded if \(m>\left(\|D\|_{\mathrm{F}}\|W\|_{\mathrm{F}}\right)^{2}\), and therefore \(m=\lceil(\|D\|_{\mathrm{F}}\|W\|_{\mathrm{F}})^{2}\rceil\), is sufficient to ensure that the covariance matrix of the output distribution \(\Psi_{m}^{\mathrm{ONN}}(\cdot,w)\stackrel{{(\mathsf{d})}}{{=}} \operatorname{Normal}(\Psi^{\mathrm{NN}}(\cdot,w),\Sigma_{\mathrm{ONN},m}^{(L)})\) in Theorem 2 is bounded in linear \(\mathsf{N}\)Ns. This is derived by using submultiplicativity of the norm (see (A.8)) and is therefore possibly a lose bound. We use the exact relation given by Corollary 2 for the covariance matrix in Theorem 2 to investigate the lowest values for \(m\) for which the covariance matrix starts being bounded. In Figure 6 we depict a linear \(\mathsf{NN}\) with constant width \(4\). We vary the values for \(\|D\|_{\mathrm{F}}\) and \(\|W\|_{\mathrm{F}}\). Upon close inspection we see that the lowest value for \(m\) seems to be \(g(x,y)\simeq\lceil(xy)^{2}/\|I_{d}\|_{\mathrm{F}}^{4}\rceil\), where \(I_{d}\) is the identity matrix of dimension \(d\), see Figure 6. Because \(\|I_{d}\|_{\mathrm{F}}=\sqrt{d}\), the value for \(m\) found numerically is \(m\simeq\lceil(\|D\|_{\mathrm{F}}\|W\|_{\mathrm{F}}/d)^{2}\rceil\).

## 6 Discussion & Conclusion

Design A, introduced in Section 3, guarantees an approximation property (Theorem 1). This is achieved through technical machinery to control the noise, even though there are nonlinear activation functions involved. This method is powerful enough to yield the universal approximation property, as NNs can be approximated arbitrarily well with \(\mathsf{ONN}\)s that

\begin{table}
\begin{tabular}{c c c c c} \hline \hline \(\#\) of additional layers & 1\&2 & 3\&4 & 5\&6 & 7\&8 \\ \hline
1 & 1 & 0 & 0 & 0 \\
2 & 1 & 1 & 0 & 0 \\
3 & 1 & 1 & 1 & 0 \\
4 & 1 & 1 & 1 & 1 \\
5 & 2 & 1 & 1 & 1 \\
6 & 2 & 2 & 1 & 1 \\ \hline \hline \end{tabular}
\end{table}
Table 1: Insertion pattern.

Figure 5: Accuracy of LeNet \(\mathsf{ONN}\)s, depending on the amount of inserted identity layers and the variance level of the \(\mathsf{ONN}\), for (a) a network with \(\tanh\) activation function and one copy, (b) a network with \(\operatorname{ReLU}\) activation function and one copy, (c) a network with linear activation function and one copy, (d) a network with \(\tanh\) activation function and two copies, (e) a network with \(\operatorname{ReLU}\) activation function and two copies, (f) a network with linear activation function and two copies.
are constructed through the first design, and NNs themselves can approximate any continuous function arbitrarily well [43, Theorem 1]. Our mathematical guarantee however, only states a sufficient number of copies required, and this number grows exponentially as the number of layers increases.

We then introduced Design B in Section 4, in which the growth of number of copies is much more benign. However, the analysis of Design B was restricted to linear NNs, and Design B might therefore not be expressive enough to have the universal approximation property. Linear NNs, or NNs with algebraic polynomials as activation functions for that matter, namely do not posses the universal approximation property. The assumption of linear activation functions did allow us to characterize the distribution of the output exactly on the flipside (Theorem 2).

In short, in this paper, we have discussed the noise present in ONNs and described a mathematical model for the noise. We also investigate the numerical implications of the the mathematical model, with a specific focus on the effects of depth (Figure 5). The proposed noise reduction schemes yield greater accuracy and the theoretical results (Theorem 1 and Corollary 3) guarantee that ONNs work just as noiseless NNs in the many copies limit. With the designs and findings of Sections 3 to 4 we have a framework to exploit known NN wisdom, as no new training is required. Further research should address optimization algorithms that take the noise of ONNs into account to investigate the regularization, generalization and minimization properties of trained ONNs.

## Acknowledgments

This research was supported by the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement no. 945045, and by the NWO Gravitation project NETWORKS under grant no. 024.002.003.

We would finally like to thank Bin Shi for advise on the noise level parameters of ONNs for our simulations. Furthermore we want to thank Albert Senen-Cerda, Sanne van Kempen, and Alexander Van Werde for feedback to a draft version of this document.
