Large-Scale Optical Neural Networks based on Photoelectric Multiplication

TL;DR

This paper introduces a scalable, high-speed, low-energy optical neural network accelerator using photoelectric multiplication, capable of reprogramming and training on the fly, with potential to surpass digital thermodynamic limits.

Contribution

It presents a novel photonic accelerator architecture based on coherent detection for large neural networks, enabling reprogrammability, training, and ultra-low-energy operation.

Findings

Achieves sub-aJ energy per MAC at GHz speeds

Reprogrammable and trainable optical neural network design

Potential to surpass thermodynamic limits for digital computation

Abstract

Recent success in deep neural networks has generated strong interest in hardware accelerators to improve speed and energy consumption. This paper presents a new type of photonic accelerator based on coherent detection that is scalable to large ($N \gtrsim 10^6$) networks and can be operated at high (GHz) speeds and very low (sub-aJ) energies per multiply-and-accumulate (MAC), using the massive spatial multiplexing enabled by standard free-space optical components. In contrast to previous approaches, both weights and inputs are optically encoded so that the network can be reprogrammed and trained on the fly. Simulations of the network using models for digit- and image-classification reveal a "standard quantum limit" for optical neural networks, set by photodetector shot noise. This bound, which can be as low as 50 zJ/MAC, suggests performance below the thermodynamic (Landauer) limit for digital irreversible computation is theoretically possible in this device. The proposed accelerator can implement both fully-connected and convolutional networks. We also present a scheme for back-propagation and training that can be performed in the same hardware. This architecture will enable a new class of ultra-low-energy processors for deep learning.

Figures

Tables

Table 1. Table 1: Layers in AlexNet [ 1 ] . Values of c in , c out subscript 𝑐 in subscript 𝑐 out c_{\rm in},c_{\rm out} are calculated from Eq. ( 7 ). Max-pooling layers after CONV1, CONV2, and CONV5 are used to reduce the image size, but the relative computational cost for these layers is negligible.

Layer	Input	Output	Kernel	Stride	MACs	$c_{\mathrm{in}}$	$c_{\mathrm{out}}$
CONV1	$227\times 227\times 3$	$55\times 55\times 96$	$11\times 11\times 96\times 3$	4	105M	93	363
(pool)	$55\times 55\times 96$	$27\times 27\times 96$	–	2	–	–	–
CONV2	$27\times 27\times 96$	$27\times 27\times 256$	$5\times 5\times 256\times 96$	1	448M	189	2400
(pool)	$27\times 27\times 256$	$13\times 13\times 256$	–	2	–	–	–
CONV3	$13\times 13\times 256$	$13\times 13\times 384$	$3\times 3\times 384\times 256$	1	150M	117	2304
CONV4	$13\times 13\times 384$	$13\times 13\times 384$	$3\times 3\times 384\times 384$	1	224M	117	3456
CONV5	$13\times 13\times 384$	$13\times 13\times 256$	$3\times 3\times 256\times 384$	1	150M	102	3456
(pool)	$13\times 13\times 256$	$6\times 6\times 256$	–	2	–	–	–
FC1	$6\times 6\times 256$	$4096$	–	–	38M	–	–
FC2	$4096$	$4096$	–	–	17M	–	–
FC2	$4096$	$1000$	–	–	4M	–	–
Total CONV layers					1.08G	132	1656
Total FC layers					59M	–	–

Table 2. Table S1: Calculation of C 0 subscript 𝐶 0 C_{0} for neural networks introduced in Sec. II.

Model	SQL $n_{\mathrm{mac}}$	Layer	$N$	$N{n}_{\mathrm{mac}}$	$C_0=N{n}_{\mathrm{mac}}{e}^2/2kT$
		1	768	3840–7680	12–24 fF
Small NN	5–10	2	100	500–1000	1.6–3.2 fF
		3	100	500–1000	1.6–3.2 fF
		1	768	384–768	1.2–2.4 fF
Large NN	0.5–1	2	1000	500–1000	1.6–3.2 fF
		3	1000	500–1000	1.6–3.2 fF

Table 3. Table S2: Gate count and transistor count for commonly used integer multipliers.

	Gate Count				Transistor Count
	8 bit	16 bit	32 bit	64 bit	8 bit	16 bit	32 bit	64 bit
Wallace and Booth	33	221	1077	4709	168	1080	5208	22680
Serial Parallel	384	1536	6144	24576	1920	7680	30720	122880
Braun Multiplier	344	1456	5984	24256	1728	7296	29952	121344
Ripple Carry Adder Based	96	384	1536	6144	480	1920	7680	30720
Vedic Multiplier	49	281	1321	5705	252	1428	6660	28644

Equations

$$x_{i}^{(k+1)}=f\Bigl{(}\sum_{j}A_{ij}^{(k)}x_{j}^{(k)}\Bigr{)}$$

$$Q_{i}=\frac{2\eta e}{\hbar\omega}\int{\mbox{Re}\bigl{[}E^{\rm(in)}(t)^{*}E_{i}^{\rm(wt)}(t)\bigr{]}{\rm d}t}\propto\sum_{j}A_{ij}x_{j}$$

$$x_{i}^{(k+1)}=f\biggl{(}\sum_{j}A_{ij}^{(k)}x_{j}^{(k)}+w_{i}^{(k)}\frac{\lVert{A^{(k)}}\rVert\lVert{x^{(k)}}\rVert}{\sqrt{N^{2}N^{\prime}}}\frac{\sqrt{N}}{\sqrt{n_{\rm mac}}}\biggr{)}$$

$$E_{\mathrm{tot}}$$

$$E_{\mathrm{mac}}$$

$${\nabla }_{A}L=({\nabla }_{Y}L)X^T,{\nabla }_{X}L=A^T({\nabla }_{Y}L)$$

$$y_{ij;k}=\sum _{i^{\prime }{j}^{\prime },l}{K}_{i^{\prime }{j}^{\prime },kl}{x}_{(s_{x}i+i^{\prime })(s_{y}j+j^{\prime });l}$$

$$E_{\mathrm{mac}}=\underset{1/c_{\mathrm{in}}}{\underbrace{\left(\frac{1}{C^{\prime }}+\frac{1}{W^{\prime }{H}^{\prime }}\right)}}{E}_{\mathrm{in}}+\underset{1/c_{\mathrm{out}}}{\underbrace{\frac{1}{K_x{K}_{y}C}}}{E}_{\mathrm{out}}$$

$$\frac{Q}{e}\sim \mathrm{Poisson}(|u{|}^2)$$

$$\frac{Q_{i}^{(\pm)}}{e}=\sum_{j}\tfrac{1}{2}(\bar{A}_{ij}\pm\bar{x}_{j})^{2}+w_{i}^{(\pm)}\Bigl{(}\sum_{j}\tfrac{1}{2}(\bar{A}_{ij}\pm\bar{x}_{j})^{2}\Bigr{)}^{1/2}$$

$${\bar{y}}_i$$

$${\bar{x}}_i^{\prime }=\alpha \text{mbox}ReLU({\bar{y}}_i^{(m)})$$

$$\bar{x}={\xi }_{x}x,{\bar{x}}^{\prime }={\xi }_x{x}^{\prime },\bar{A}={\xi }_{A}A$$

$$x^{\prime}_{i}=2\alpha\xi_{A}\,\mbox{ReLU}\biggl{[}A_{ij}x_{j}+\frac{w_{i}}{2}\sqrt{\frac{\lVert{A_{i}}\rVert^{2}}{\xi_{x}^{2}}+\frac{\lVert{x}\rVert^{2}}{\xi_{A}^{2}}}\biggr{]}$$

$$n_{\mathrm{tot}}^{(x)}=N^{\prime }\Vert x{\Vert }^2{\xi }_x^2,n_{\mathrm{tot}}^{(A)}=\Vert A{\Vert }^2{\xi }_A^2$$

$$n_{\mathrm{mac}}^{(x)}=\frac{\Vert x{\Vert }^2}{N}{\xi }_x^2,n_{\mathrm{mac}}^{(A)}=\frac{\Vert A{\Vert }^2}{N{N}^{\prime }}{\xi }_A^2$$

$$x^{\prime}_{i}=\mbox{ReLU}\biggl{[}A_{ij}x_{j}+w_{i}\frac{\lVert{A}\rVert\lVert{x}\rVert}{2\sqrt{NN^{\prime}}}\sqrt{\frac{1}{n_{\rm mac}^{(x)}}+\frac{1}{n_{\rm mac}^{(A)}}}\biggr{]}$$

$$x^{\prime}_{i}=\mbox{ReLU}\biggl{[}A_{ij}x_{j}+w_{i}\frac{\lVert{A}\rVert\lVert{x}\rVert}{\sqrt{N^{2}N^{\prime}}}\frac{\sqrt{N}}{\sqrt{n_{\rm mac}}}\biggr{]}$$

$${\bar{C}}_{ij}=2\sum _l{\bar{A}}_{il}{\bar{B}}_{lj}+w_{ij}\sqrt{\Vert {\bar{A}}_{i,:}{\Vert }^2+\Vert {\bar{B}}_{:,j}{\Vert }^2}$$

$$C_{ij}=\sum _l{A}_{il}{B}_{lj}+\frac{w_{ij}}{2}\sqrt{\frac{\Vert A_{i,:}{\Vert }^2}{{\xi }_A^2}+\frac{\Vert B_{:,j}{\Vert }^2}{{\xi }_B^2}}$$

$$C_{ij}=\sum _l{A}_{il}{B}_{lj}+\frac{w_{ij}}{2}\frac{\Vert A\Vert \Vert B\Vert }{\sqrt{mnk}}\sqrt{\frac{1}{n_{\mathrm{mac}}^{(A)}}+\frac{1}{n_{\mathrm{mac}}^{(B)}}}$$

$$C_{ij}=\sum _l{A}_{il}{B}_{lj}+w_{ij}\frac{\Vert A\Vert \Vert B\Vert }{\sqrt{mnk}}\frac{1}{\sqrt{n_{\mathrm{mac}}}}$$

$${\bar{y}}_i=2{\bar{A}}_{ij}{\bar{x}}_j+w_i\sqrt{\Vert {\bar{A}}_i{\Vert }^2+\Vert \bar{x}{\Vert }^2+2⟨\Delta n_e^2⟩}$$

$$x^{\prime}_{i}=\mbox{ReLU}\biggl{[}A_{ij}x_{j}+w_{i}\frac{\lVert{A}\rVert\lVert{x}\rVert}{2\sqrt{NN^{\prime}}}\sqrt{\frac{1}{n_{\rm mac}^{(x)}}+\frac{1}{n_{\rm mac}^{(A)}}+\frac{2\langle\Delta n_{e}^{2}\rangle}{Nn_{\rm mac}^{(A)}n_{\rm mac}^{(x)}}}\biggr{]}$$

$$x^{\prime}_{i}=\mbox{ReLU}\biggl{[}A_{ij}x_{j}+w_{i}\frac{\lVert{A}\rVert\lVert{x}\rVert}{\sqrt{NN^{\prime}}}\frac{1}{\sqrt{n_{\rm mac}}}\sqrt{1+\frac{2\langle\Delta n_{e}^{2}\rangle}{Nn_{\rm mac}}}\biggr{]}$$

$$\frac{⟨\Delta n^2{⟩}_{\mathrm{th}}}{⟨\Delta n^2{⟩}_{\mathrm{shot}}}=\frac{2⟨\Delta n_e^2⟩}{N{n}_{\mathrm{mac}}}=\frac{C}{N{n}_{\mathrm{mac}}{e}^2/2kT}\equiv \frac{C}{C_0}$$

Full text

Research Laboratory of Electronics, MIT, 50 Vassar Street, Cambridge, MA 02139, USA

In recent years, deep neural networks have tackled a wide range of problems including image analysis [1], natural language processing [2], game playing [3], physical chemistry [4], and medicine [5]. This is not a new field, however. The theoretical tools underpinning deep learning have been around for several decades [6, 7, 8]; the recent resurgence is driven primarily by (1) the availability of large training datasets [9], and (2) substantial growth in computing power [10] and the ability to train networks on GPUs [11]. Moving to more complex problems and higher network accuracies requires larger and deeper neural networks, which in turn require even more computing power [12]. This motivates the development of special-purpose hardware optimized to perform neural-network inference and training [13].

To outperform a GPU, a neural-network accelerator must significantly lower the energy consumption, since the performance of modern microprocessors is limited by on-chip power [14]. In addition, the system must be fast, programmable, scalable to many neurons, compact, and ideally compatible with training as well as inference. Application-specific integrated circuits (ASICs) are one obvious candidate for this task. State-of-the-art ASICs can reduce the energy per multiply-and-accumulate (MAC) from 20 pJ/MAC for modern GPUs [15] to around 1 pJ/MAC [16, 17]. However, ASICs are based on CMOS technology and therefore suffer from the interconnect problem—even in highly optimized architectures where data is stored in register files close to the logic units, a majority of the energy consumption comes from data movement, not logic [13, 16]. Analog crossbar arrays based on CMOS gates [18] or memristors [19, 20] promise better performance, but as analog electronic devices, they suffer from calibration issues and limited accuracy [21].

Photonic approaches can greatly reduce both the logic and data-movement energy by performing (the linear part of) each neural-network layer in a passive, linear optical circuit. This allows the linear step is performed at high speed with no energy consumption beyond transmitter and receiver energies. Optical neural networks based on free-space diffraction [22] have been reported, but require spatial light modulators or 3D-printed diffractive elements, and are therefore not rapidly programmable. Nanophotonic circuits are a promising alternative [23, 24], but the footprint of directional couplers and phase modulators makes scaling to large ($N\geq 1000$) numbers of neurons very challenging. To date, the goal of a large-scale, rapidly reprogrammable photonic neural network remains unrealized.

This paper presents a new architecture based on coherent (homodyne) detection that is fast, low-power, compact, and readily scalable to large ($N\gtrsim 10^{6}$) numbers of neurons. In contrast to previous schemes, here we encode both the inputs and weights in optical signals, allowing the weights to be changed on the fly at high speed. Synaptic connections (matrix-vector products) are realized by the quantum photoelectric multiplication process in the homodyne detectors. Our system is naturally adapted to free-space optics and can therefore take advantage of the massive spatial multiplexing possible in free-space systems [25, 26] and the high pixel density of modern focal-plane arrays [27] to scale to far more neurons than can be supported in nanophotonics or electronic cross-bar arrays. The optical energy consumption is subject to a fundamental standard quantum limit (SQL) arising from the effects of shot noise in photodetectors, which lead to classification errors. Simulations based on MNIST neural networks empirically show the SQL can be as low as 50–100 zJ/MAC. Using realistic laser, modulator, and detector energies, performance at the sub-fJ/MAC level should be possible with present technology. The optical system can be used for both fully-connected and convolutional layers. Finally, backpropagation is straightforward to implement in our system, allowing both inference and training to be performed in the same optical device.

Coherent Matrix Multiplier

Fig. 1 illustrates the device. A deep neural network is a sequence of $K$ layers (Fig. 1(a)), each consisting of a matrix multiplication $\vec{x}\rightarrow A\vec{x}$ (synaptic connections) and an element-wise nonlinearity $x_{i}\rightarrow f(x_{i})$ (activation function); thus the input into the $(k+1)^{\rm th}$ layer is related to the $k^{\rm th}$ layer input by:

[TABLE]

For a given layer, let $N$ and $N^{\prime}$ be the number of input and output neurons, respectively. Input (output) data are encoded temporally as $N$ ($N^{\prime}$) pulses on a single channel as shown in Fig. 1(b). This encoding, reminiscent of the Coherent Ising Machine [28, 29, 30], contrasts with other approaches used for neural networks, which encode inputs in separate spatial channels [23, 24, 22]. As there are $NN^{\prime}$ weights for an $N^{\prime}\times N$ fully-connected matrix, the weights enter on $N^{\prime}$ separate channels, each carrying a single matrix row encoded in time. Input data is optically fanned out to all $N^{\prime}$ channels, and each detector functions as a quantum photoelectric multiplier, calculating the homodyne product between the two signals (inset). As long as both signals are driven from the same coherent source and the path-length difference is less than the coherence length, the charge $Q_{i}$ accumulated by homodyne receiver $i$ is:

[TABLE]

Here $E^{\rm(in)}(t)$ and $E_{i}^{\rm(wt)}(t)$ are the input and weight fields for receiver $i$, which are taken to be sequences of pulses with amplitudes proportional to $x_{j}$ and $A_{ij}$, respectively ($x_{j},A_{ij}\in\mathbb{R}$). Thus each receiver performs a vector-vector product between $\vec{x}$ and a row $\vec{A}_{i}$ of the weight matrix; taken together, the $N^{\prime}$ electronic outputs give the matrix-vector product $A\vec{x}$. Fields are normalized so that power is given by $P(t)=|E(t)|^{2}$, and $\eta$ is the detector efficiency. A serializer reads out these values one by one, applies the nonlinear function $f(\cdot)$ in the electrical domain, and outputs the result to a modulator to produce the next layer’s inputs.

The balanced homodyne detector in Fig. 1(b) (inset) combines the advantages of optics and electronics: it can process data encoded at extremely high speeds, limited only by the bandwidth of the beamsplitter ($\gtrsim$THz) and the (optical) bandwidth of the photodetectors (typically $\gtrsim$ 100 nm, or $\gtrsim$ 10 THz). The electrical bandwidth can be much slower, since only the integrated charge is measured. Finally, the present scheme avoids the need for low-power nonlinear optics that is a major stumbling block in all-optical logic [31]: since the output is electrical, the dot product $A_{ij}x_{j}$ can be computed at extremely low power (sub-fJ/MAC) using standard non-resonant components (photodiodes) that are CMOS-compatible and scalable to arrays of millions.

Previous approaches used optoelectronics (photodiodes, lasers, amplifiers) both to sum neuron inputs [32, 24] and to generate nonlinearity or spiking dynamics [33, 34, 35, 36, 37]; here, thanks to the optical weight encoding, the synaptic weighting itself is performed opto-electronically.

Coherent detection greatly simplifies the setup compared to alternative approaches. With a given set of weight inputs, the network in Fig. 1(b) requires $N$ input pulses and $N^{\prime}$ detectors to perform a matrix-vector operation with $NN^{\prime}$ MACs, performing an operation that should scale quadratically with size using only linear resources. This is in contrast to electrical approaches that require quadratic resources ($NN^{\prime}$ floating-point operations total). The (optical) energy consumption of nanophotonic systems [23, 24] also scales linearly for the same operation; however, the circuit is much more complex, requiring $O(NN^{\prime})$ tunable phase shifters [38, 39] or ring resonators [24], which becomes very challenging to scale beyond several hundred channels and may be sensitive to propagation of fabrication errors. The main caveat to our system is the need to generate the weights in the first place, which imposes an energy cost that does scale quadratically. However, in many cases (particularly in data centers) neural networks are run simultaneously over large batches of data, so with appropriate optical fan-out, the cost of the weights can be amortized over many clients. Put another way, running the neural network on data with batch size $B$, we are performing a matrix-matrix product $Y_{N^{\prime}\times B}=A_{N^{\prime}\times N}X_{N\times B}$, which requires $N^{\prime}NB$ MACs, with an energy cost that should scale as $O(N^{\prime}N)+O(N^{\prime}B)+O(NB)$ rather than $O(N^{\prime}NB)$.

Deep Learning at the Standard Quantum Limit

As energy consumption is a primary concern in neuromorphic and computing hardware generally [14], an optical approach must outperform electronics by a large factor to justify the investment in a new technology. In addition, optical systems must show great potential for improvement, ideally by many orders of magnitude, to allow continued scaling beyond the physical limits of Moore’s Law. Thus two questions are relevant: (1) the fundamental, physical limits to the energy consumption, and (2) the energy consumption of a practical, near-term device using existing technology.

The fundamental limit stems from quantum-limited noise. In an electrical signal, energy is quantized at a level $E_{\rm el}=h/\tau_{\rm el}$, where $\tau_{\rm el}\sim 10^{-10}\,{\rm s}$ is the signal duration. Optical energy is quantized at a level $E_{\rm opt}=h/\tau_{\rm opt}$, where $\tau_{\rm opt}\equiv c/\lambda\sim(2$–$5)\times 10^{-15}\,{\rm s}$, which is $10^{4}$–$10^{5}$ times higher. As a result, $E_{\rm opt}\gg kT\gg E_{\rm el}$ and electrical signals can be treated in a classical limit governed by thermal noise, while optical signals operate in a zero-temperature quantum limit where vacuum fluctuations dominate. These fluctuations are read out on the photodetectors, where the photoelectric effect [40] produces a Poisson-distributed photocurrent [41, 42]. While the photocurrents are subtracted in homodyne detection, the fluctuations add in quadrature, and Eq. (1) is replaced by (See Supp. Sec. S3 for derivation and assumptions):

[TABLE]

Here the $w_{i}^{(k)}\sim N(0,1)$ are Gaussian random variables, $\lVert{\cdot}\rVert$ is the $L_{2}$ norm, and $n_{\rm mac}$ is the number of photons per MAC, related to the total energy consumption of the layer by $n_{\rm tot}=NN^{\prime}n_{\rm mac}$.

The noise term in Eq. (3) scales as $n_{\rm mac}^{-1/2}$, and therefore the signal-to-noise ratio (SNR) of each layer will scale as ${\rm SNR}\propto n_{\rm mac}$. Since noise adversely affects the network’s performance, one expects that the energy minimum should correspond to the value of $n_{\rm mac}$ at which the noise becomes significant. To quantify this statement, we perform benchmark simulations using a collection of neural networks trained on the MNIST (digit recognition) dataset. While MNIST digit classification is a relatively easy task [13], the intuition developed here should generalize to more challenging problems. Data for two simple networks are shown in Fig. 2, both having a 3-layer, fully-connected topology (Fig. 2(a)). In the absence of noise, the networks classify images with high accuracy, as the example illustrates (Fig. 2(b)).

As Fig. 2(c) shows, the error rate is a monotonically decreasing function of $n_{\rm mac}$. The two asymptotic limits correspond to the noiseless case ($n_{\rm mac}\rightarrow\infty$, which returns the network’s canonical accuracy), and the noise-dominated case ($n_{\rm mac}\rightarrow 0$, where the network is making a random guess). Of interest to us is the cutoff point, loosely defined as the lowest possible energy at which the network returns close to its canonical accuracy (for example within a factor of $2\times$, see dashed lines in Fig. 2(c)). This is around 0.5–1 aJ (5–10 photons) for the small network (inner layer size $N=100$), and 50–100 zJ (0.5–1 photon) for the large network (inner layer size $N=1000$). (Note that this is per MAC, the number of photons per detector $Nn_{\rm mac}$ is typically $\gg 1$.) This bound stems from the standard quantum limit (SQL): the intrinsic uncertainty of quadrature measurements on coherent states [43], which is temperature- and device-independent. This should be viewed as an absolute lower bound for the energy consumption of neural networks of this type; although the use of squeezed light allows one to reach sensitivity below the SQL [44, 45], this requires squeezing all inputs (including vacuum inputs in optical fan-out) which will likely lead to a net increase in overall energy consumption (squeezing injects an average of $\sinh^{2}(\eta)$ photons per pulse, where $\eta$ is the squeezing parameter [41], which will substantially increase $n_{\rm mac}$).

The SQL is network-dependent, and not all layers contribute equally. For each MAC, we have ${\rm SNR}\propto n_{\rm mac}$; however, the signal adds linearly while the errors add in quadrature. As a result, the larger network is more resilient to individual errors because each output is averaging over more neurons. Moreover, the solid curves in Fig. 2(c) are restricted to the case when $n_{\rm mac}$ is the same for all layers. The dashed lines show the error rate in a fictitious device where quantum-limited noise is only present in a particular layer. For the large network, a smaller $n_{\rm mac}$ can be tolerated in the second layer, suggesting that better performance could be achieved by independently tuning the energy for each layer. Moreover, just as neural networks can be “co-designed” to achieve high accuracy on limited bit-precision hardware [13], changes to the training procedure (e.g. injecting noise to inner layers, a technique used to reduce generalization error [46, 47]) may further improve performance at low powers.

Quantum limits to computational energy efficiency in photonics are not unique to neural networks. In digital photonic circuits based on optical bistability [48], vacuum fluctuations lead to spontaneous switching events that limit memory lifetime and gate accuracy [49, 50]. However, these effects require bistability at the attojoule scale [49, 51], which is well out of the reach of integrated photonics (although recent developments are promising [52, 53, 54]). By contrast, neural networks are analog systems so the quantum fluctuations set a meaningful limit on efficiency even though no attojoule-scale optical nonlinearities are employed.

Energy Budget

Viewing the neural network as an analog system with quantum-limited performance shifts the paradigm for comparing neural networks. Fig. 3(a) shows the standard approach: a scatterplot comparing error rate with number of MACs, a rough proxy for time or energy consumption [12, 13]. There is a tradeoff between size and accuracy, with larger networks requiring more operations but also giving better accuracy. In the SQL picture, each point becomes a curve because now we are free to vary the number of photons per MAC, and the energy bound is set by the total number of photons, not the number of MACs. Fig. 3(b) plots the error rate as a function of photon number for the networks above. While the general tradeoff between energy and accuracy is preserved, there are a number of counterintuitive results. For example, according to Fig. 3(a), networks 1 and 2 have similar performance but the first requires $8\times$ more MACs, so under a conventional analysis, network 2 would always be preferred. However, Fig. 3(b) indicates that network 1 has better performance at all energy levels. This is because network 1 is less sensitive to shot noise due to averaging over many neurons, and therefore can be operated at lower energies, compensating for the increased neuron count. The same apparent paradox is seen with networks 3 and 4. This suggests that, in a quantum-limited scenario, reducing total energy may not be as simple as reducing the number of operations.

The total energy budget depends on many factors besides the SQL. Fig. 3(c) plots energy per MAC as a function of the average number of input neurons per layer $N$, a rough “size” of the neural network. The SQL data are plotted for the eight networks in Fig. 3(a-b), and the corresponding dashed line is an empirical fit. Note that the SQL is an absolute lower bound, assumes perfect detectors, and only counts input optical energy. In a realistic device, this curve is shifted up by a factor $(\eta_{d}\eta_{c}\eta_{s}\beta_{\rm mod})^{-1}$, where $\eta_{d}$, $\eta_{c}$, and $\eta_{s}$ are the detector, coupling, and source (laser) efficiencies and $\beta_{\rm mod}$ is the modulator launch efficiency [55]; these are all close enough to unity in integrated systems [26, 56, 57, 58] that the factor is $\lesssim 10$.

Another key factor is the detector electronics. The homodyne signal from each neuron needs to be sent through a nonlinear function $y_{i}\rightarrow f(y_{i})$ and converted to the optical domain using a modulator (Fig. 1(b)). The most obvious way to do this is to amplify and digitize the signal, perform the function $f(\cdot)$ in digital logic, serialize the outputs, convert back to analog, and send the analog signal into the modulator. Transimpedance amplifiers designed for optical interconnects operate at the $\sim 100$ fJ range [59, 26], while ADCs in the few-pJ/sample regime are available [60] and simple arithmetic (for the activation function) can be performed at the pJ scale [15, 16, 17]. Modulators in this energy range are standard [56, 57, 59]. Thus a reasonable near-term estimate would be few-pJ/neuron; this figure is divided by the number of inputs per neuron to give the energy per MAC (solid green curve in Fig. 3(c)). This few-pJ/neuron figure includes both optical and electrical energy: even though only a fraction of the energy is optical, the optical signal will be large compared to both shot noise (Eq. (3)) and amplifier Johnson noise $\langle\Delta n_{e}\rangle_{\rm rms}\sim 10^{3}$ [61], so noise will not significantly degrade the network’s performance.

A much more aggressive goal is 1 fJ/neuron (dashed green curve). This figure is out of reach with current technology, but research into fJ/bit on-chip interconnects may enable it in the future [61, 26]. A range of modulator designs support few-fJ/bit operation [62, 63, 64, 65]. On-chip interconnects also require photodetectors with ultra-low (fF) capacitance, so that a femtojoule of light produces a detectable signal without amplification [61, 26]; such detectors have been realized with photonic crystals [66], plasmon antennas [67, 68], and nanowires [69]. By eliminating the amplifier, ultrasmall “receiverless” detectors avoid its $\sim 100$ fJ energy cost as well as the Johnson noise associated with the amplifier. (Johnson noise still leads to fluctuations in the capacitor charge (kTC noise) that go as $\langle\Delta n_{e}\rangle_{\rm rms}=\sqrt{kTC}/e\approx 12\sqrt{C/{\rm fF}}$ [70], but for small detectors shot noise will dominate, see Supp. Sec. S4). Since 1 fJ/neuron is below the energy figures for ADCs, it would require well-designed analog electronics (for the nonlinear activation function) and very tight integration between detector, logic, and modulator [26]. At these energies, shot noise is also non-negligible and the SQL becomes relevant, but as mentioned above, due to optical inefficiencies the SQL will likely be relevant at higher energies as well.

For context, the $\sim$1 pJ/MAC figure [15, 16, 17] for state-of-the-art ASICs is shown in Fig. 3(c). Energy consumption in non-reversible logic gates is bounded by the Landauer (thermodynamic) limit $E_{\rm op}=kT\log(2)\approx 3\,{\rm zJ}$ [71]. While multiply-and-accumulate is technically a reversible operation, all realistic computers implement it using non-reversible binary gates, so Landauer’s principle applies. A 32-bit multiplication [72, 73] requires approximately $10^{3}$ binary gates (see Supp. Sec. S4) and each bit operation consumes at least $kT\log(2)$, giving a limit $E_{\rm mac}\geq 3\,{\rm aJ}$ (dotted line in Fig. 3(c)). This is already higher than the SQL for the larger networks with $N\geq 100$. The optical neural network can achieve sub-Landauer performance because (1) it operates in analog, avoiding the overhead of many bit operations per multiplication, and (2) the matrix product is performed through optical interference, which is reversible and not subject to the bound. To understand the second point, recall that homodyne detection computes the dot product via the polarization identity: $\vec{u}\cdot\vec{v}=\tfrac{1}{4}(\lVert\vec{u}+\vec{v}\rVert^{2}-\lVert\vec{u}-\vec{v}\rVert^{2})$. Optical interference, the reversible element that breaks Landauer’s assumption, is needed to convert the signals representing $\vec{u}$ and $\vec{v}$ to $\vec{u}\pm\vec{v}$ before squaring on the detectors and subtracting.

A final consideration is the electrical energy required to generate the weights. There is one weight pulse per MAC, so at the minimum this will be 1 fJ/MAC for the modulator, and may rise above 1 pJ/MAC once the driver electronics and memory access are included. However, once the optical signal is generated, it can be fanned out to many neural networks in parallel, reducing this cost by a factor of $B$, the batch size. Large batch sizes should enable this contribution to $E_{\rm mac}$ to reach the few-femtojoule regime, and potentially much lower.

Training and Convolutions with Optical GEMM

As discussed previously, the optical unit in Fig. 1(b) performs a matrix-vector product, and running multiple units in parallel with the same set of weights performs a general matrix-matrix product (GEMM), a key function in the Basic Linear Algebra Subprograms (BLAS) [74]. Fig. 4(a-b) shows a schematic for an optical GEMM unit based on homodyne detection inspired by the neural-network concept. The inputs are two matrices $(M_{1})_{m\times k}$ and $(M_{2})_{n\times k}$, encoded into optical signals on the 1D red (blue) integrated photonic transmitter arrays. Cylindrical lenses map these inputs to rows (columns) of the 2D detector array. From the accumulated charge at each pixel, one can extract the matrix elements of the product $(M_{1}M_{2}^{T})_{m\times n}$. This operation requires $m\cdot n\cdot k$ MACs, and the total energy consumption (and energy per MAC) are:

[TABLE]

where $E_{\rm in}$, $E_{\rm out}$ are the transmitter and receiver energy requirements, per symbol, which include all optical energy plus electronic driving, serialization, DAC/ADC, etc. If all matrix dimensions $(m,n,k)$ are large, significant energy savings per MAC are possible if $E_{\rm in},E_{\rm out}$ can be kept reasonably small.

We saw above that the optical system could be used for neural-network inference. When running a batch of $B$ instances $X=[x_{1}\ldots x_{B}]$, the output $Y=[y_{1}\ldots y_{B}]$ can be computed through the matrix-matrix product $Y=AX$. In fully-connected layers, training and back-propagation also rely heavily on GEMM. The goal of training is to find the set of weights $A^{(k)}$ that minimize the loss function $L$, which characterizes the inaccuracy of the model. Training typically proceeds by gradient-based methods. Since the loss depends on the network output, we start at the final layer and work backward, a process called back-propagation [7, 8]. At each layer, we compute the gradient $(\nabla_{A}L)_{ij}=\partial L/\partial A_{ij}$ from the quantity $(\nabla_{Y}L)_{ij}=\partial L/\partial Y_{ij}$, and propagate the derivative back to the input $(\nabla_{X}L)_{ij}=\partial L/\partial X_{ij}$ (Fig. 4(c)). These derivatives are computed from the chain rule and can be written as matrix-matrix multiplications:

[TABLE]

Once the derivative has been propagated to $\nabla_{X^{(k)}}L$ (for layer $k$) we use the chain rule to compute $\nabla_{Y^{(k-1)}}L=f^{\prime}(\nabla_{X^{(k)}}L)$ and proceed to the previous layer. In this way, we sequentially compute the derivatives $\nabla_{A^{(k)}}L$ at each layer in the neural network.

In addition to fully-connected layers, it is also possible to run convolutional layers on the optical GEMM unit by employing a “patching” technique [75]. In a convolutional layer, the input $x_{ij;k}$ is a $W\times H$ image with $C$ channels. This is convolved to produce an output $y_{ij;k}$ of dimension $W^{\prime}\times H^{\prime}$ with $C^{\prime}$ channels [13]:

[TABLE]

Here $K_{i^{\prime}j^{\prime},kl}$ is the convolution kernel, a 4-dimensional tensor of size $K_{x}\times K_{y}\times C^{\prime}\times C$, and $(s_{x},s_{y})$ are the strides of the convolution. Naïvely vectorizing Eq. (6) and running it as a fully-connected matrix-vector multiply is very inefficient because the resulting matrix is sparse and contains many redundant entries. Patching expresses the image as a matrix $X$ of size $K_{x}K_{y}C\times W^{\prime}H^{\prime}$, where each column corresponds to a vectorized $K_{x}\times K_{y}$ patch of the image (Fig. 4(d)). The elements of the kernel are rearranged to form a (dense) matrix $K$ of size $C^{\prime}\times K_{x}K_{y}C$. Eq. (6) can then be computed by taking the matrix-matrix product $Y=KX$, which has size $C^{\prime}\times W^{\prime}H^{\prime}$. On virtually any microprocessor, GEMM is a highly optimized function with very regular patterns of memory access; the benefits of rewriting the convolution as a GEMM greatly outweigh the redundancy of data storage arising from overlapping patches [75]. The time required to rearrange the image as a patch matrix is typically very small compared to the time to compute the GEMM [76] (and can be further reduced if necessary with network-on-chip architectures [77] or optical buffering [78]); therefore, by accelerating the GEMM, the optical matrix multiplier will significantly improve the speed and energy efficiency of convolutional layers. Note also that, since we are performing the convolution as a matrix-matrix (rather than matrix-vector) operation, it is possible to obtain energy savings even without running the neural network on large batches of data. Computing the convolution requires $W^{\prime}H^{\prime}K_{x}K_{y}C^{\prime}C$ MACs. Following Eq. (4), the energy per MAC (not including memory rearrangement for patching) is:

[TABLE]

The coefficients $c_{\rm in}=(1/C+1/W^{\prime}H^{\prime})^{-1}$ and $c_{\rm out}=K_{x}K_{y}C$ govern the energy efficiency when we are limited by input / output energies (transmitter / receiver and associated electronics). Since reading a 32-bit register takes $\sim$pJ of energy [13], a reasonable lower bound for near-term systems is $E_{\rm in},E_{\rm out}\gtrsim{\rm pJ}$. Thus it is essential that $c_{\rm in},c_{\rm out}\gg 1$ for the energy performance of the optical system to beat an ASIC ($\sim$pJ/MAC).

As a benchmark problem, we consider AlexNet [1], the first convolutional neural network to perform competitively at the ImageNet Large-Scale Visual Recognition Challenge [9]. AlexNet consists of 5 convolutional (CONV) layers and 3 fully-connected (FC) layers, and consistent with deep neural networks generally, the majority of the energy consumption comes from the CONV layers [13]. Table 1 gives the layer dimensions and the values of $c_{\rm in},c_{\rm out}$ for the CONV layers in AlexNet [1]. The MAC-weighted averages for all layers are $\langle c_{\rm in}\rangle>100$ and $\langle c_{\rm out}\rangle>1000$. Thus, even under extremely conservative assumptions of $E_{\rm in},E_{\rm out}\gtrsim 100\,{\rm pJ}$ (comparable to DRAM read energies [13, 14]), it is still possible to achieve sub-pJ/MAC performance.

More advanced technology, such as few-fJ optical interconnects [26], may significantly reduce $E_{\rm in}$ and $E_{\rm out}$, and therefore the energy per MAC. However, the performance is still fundamentally limited by detector shot noise (e.g. Eq. (3) for FC layers). Supp. Sec. S3 extends the shot-noise analysis to the case of matrix-matrix products needed for the convolutional case. Using a pre-trained AlexNet model (see Methods for details), Fig. 5(b) shows the top-10 accuracy on the ImageNet validation set as a function of the number of photons per MAC $n_{\rm mac}$. Consistent with Fig. 2(c), there are two limits: $n_{\rm mac}\ll 1$ corresponds to the random guess regime with 99% error rate (for top-10 accuracy with 1,000 classes), while $n_{\rm mac}\gg 1$ recovers the accuracy of the noiseless model.

The dashed lines in Fig. 5(b) show the fictitious case where noise is present in only a single layer, while the solid green line corresponds to the case where all layers have noise and $n_{\rm mac}$ is the same for each layer. Not all layers contribute equally to the noise: CONV1 is the most sensitive, requiring $n_{\rm mac}\gtrsim 20$, while the deeper layers (particularly the fully-connected layers) can tolerate much lower energies $n_{\rm mac}\gtrsim 1$. Since the SNR is related to the total power received, which scales as $c_{\rm out}n_{\rm mac}$ for the convolutional layers ($c_{\rm out}$ pulses per detector), it is not surprising that the deeper layers, which have a larger $c_{\rm out}$, are less sensitive to quantum noise. The SQL obtained for AlexNet ($n_{\rm mac}\gtrsim 20$ or $E_{\rm mac}\gtrsim 3\,$aJ) is slightly larger than that from the MNIST networks in Fig. 2(c), but of the same order of magnitude, suggesting that the SQL is somewhat problem-dependent.

It is worth contrasting the optical GEMM to more familiar optical convolvers. It has long been known that 2D convolutions can be performed with optical Fourier transforms [79, 80, 81]. However, this technique suffers from two significant drawbacks. First, it employs spatial light modulators, which limits the speed at which the kernel can be reprogrammed. In addition, optics performs a single-channel ($C=C^{\prime}=1$) convolution, and while extending to multiple output channels is possible by tiling kernels [82], multiple input and output channels may be difficult.

In contrast to free-space and fully-integrated approaches, the optical GEMM leverages the complementary strengths of both free-space and integrated photonics. Integrated photonics is an ideal platform for realizing the transmitters, as these employ a large number of fast (GHz) modulators on chip. On-chip integration allows scaling to large arrays with control over the relative phase of each output beam (a capability exploited in recent chip-based phased arrays for beam steering [83, 84, 85]). Free-space propagation provides an essential third dimension, which enables high bandwidths at moderate clock frequencies [26] and data fan-out patterns that are difficult to implement on a 2D photonic chip. However, having a free-space element leads to issues with phase stability and aberrations. Since the transmitters are integrated, it is the relative phase between the beam paths that drifts (on timescales long compared to a computation), and this can be stabilized with a single feedback loop to the overall transmitter phase, a small constant overhead that does not scale with matrix size. To correct for geometric aberrations and minimize crosstalk between detectors, multi-lens arrangements can be used, a standard practice in high-resolution imaging systems [86]. Supp. Sec. S2 presents an optical design and analysis using Zemax® simulation software supporting the hypothesis that a $10^{3}\times 10^{3}$ optical GEMM is achievable.

Discussion

This paper has presented a new architecture for optically accelerated deep learning that is scalable to large problems and can operate at high speeds with low energy consumption. Our approach takes advantage of the photoelectric effect, via the relation $I\propto|E|^{2}$, to compute the required matrix products opto-electronically without need for an all-optical nonlinearity, a key difficulty that has hobbled conventional approaches to optical computing [31]. Since the device can be constructed with free-space optical components, it can scale to much larger sizes than purely nanophotonic implementations [23], being ultimately limited by the size of the detector array ($N\gtrsim 10^{6}$).

A key advantage to this scheme is that the multiplication itself is performed passively by optical interference, so the main speed and energy costs are associated with routing data into and out of the device. For a matrix multiplication $C_{m\times n}=A_{m\times k}B_{k\times n}$, the input/output (I/O) energy scales as $O(mk)+O(nk)+O(mn)$, while the number of MACs scales as $O(mnk)$. For moderately large problems found in convolutional neural-network layers ($m,n,k\geq 100$) with moderate I/O energies ($\sim$pJ), performance in the $\sim$10 fJ/MAC range should be feasible, which is 2–3 orders of magnitude smaller than for state-of-the-art CMOS circuits [15, 16, 17]. Advances in optical interconnects [62, 56, 57] may reduce the I/O energies by large factors [26], translating to further improvements in energy per MAC.

The fundamental limits to a technology are important to its long-term scaling. For the optical neural network presented here, detector shot noise presents a standard quantum limit to neural network energy efficiency [43]. Because this limit is physics-based, it cannot be engineered away unless non-classical states of light are employed [44, 45]. To study the SQL in neural networks, we performed Monte Carlo simulations on pre-trained models for MNIST digit recognition (fully-connected) and ImageNet image classification (convolutional). In both cases, network performance is a function of the number of photons used, which sets a lower bound on the energy per MAC. This bound is problem- and network-dependent, and for the problems tested in this paper, lies in the range $50\,{\rm zJ}$–$5\,{\rm aJ}$/MAC. By contrast, the Landauer (thermodynamic) limit for a digital processor is 3 aJ/MAC (assuming 1,000 bit operations per MAC [72, 73]); sub-Laudauer performance is possible because the multiplication is performed through optical interference, which is reversible and not bounded by Landauer’s principle.

Historically, the exponential growth in computing power has driven advances in machine learning by enabling the development of larger, deeper, and more complex models [11, 13, 16, 17]. As Moore’s Law runs out of steam, photonics may become necessary for continued growth in processing power—not only for interconnects [26], but also for logic. The architecture sketched in this paper promises significant short-term performance gains over state-of-the-art electronics, with a long-term potential, bounded by the standard quantum limit, of many orders of magnitude of improvement.

Methods

Neural-network performance was computed using Monte Carlo simulations. For fully-connected layers, Eq. (3) was used, while for convolutional layers, the convolution was performed by first forming the patch matrix (Fig. 4(d)) and performing the matrix-matrix multiplication (noise model discussed in Supp. Sec. S3). The weights for the fully-connected MNIST neural networks were trained on a GPU using TensorFlow. A pretrained TensorFlow version of AlexNet (available online at Ref. [87]) was modified to implement the quantum noise model and used for ImageNet classification. Simulations were performed on an NVIDIA Tesla K40 GPU.

Acknowledgements

R.H. is supported by an IC Postdoctoral Research Fellowship at MIT, administered by ORISE through U.S. DOE and ODNI. L.B. is supported by a Doctoral Postgraduate Scholarship from the Natural Sciences and Engineering Research Council of Canada (NSERC). D.E. and M.S. acknowledge support from the U.S. ARO through the ISN at MIT (no. W911NF-18-2-0048). The authors acknowledge John Peurifoy (MIT) for training a number of the MNIST neural-network models, and NVIDIA Corporation for the donation of the Tesla K40 GPU used in this research. We are grateful to Joel Emer (NVIDIA / MIT) and Vivienne Sze (MIT) for helpful discussions.

Correspondence and requests for materials should be addressed to R.H. (email: rhamerly@mit.edu).

[Supplementary] Large-Scale Optical Neural Networks based on Photoelectric Multiplication

S1 Homodyne Product Implementation Details

Using optical homodyne detection (Fig. S1(a)), it is possible to obtain a signal proportional to the product of electric field amplitudes originating from two coherent, in-phase optical sources at different spatial locations. The input beams of electric field amplitudes $x_{1}$ and $A_{11}$ travel through a 50:50 beamsplitter, interfere, and outputs $I_{+}$ and $I_{-}$ are detected, where $I_{+}=\tfrac{1}{2}|x_{1}+A_{11}|^{2}$ and $I_{-}=\tfrac{1}{2}|x_{1}-A_{11}|^{2}$. The difference of the photocurrents $I_{+}-I_{-}$ is proportional to the real part of the product of the incident electric field amplitudes $2\mbox{Re}[A_{11}^{*}x_{1}]$. If all field amplitudes are real, this returns the product $2A_{11}x_{1}$.

Alternatively, to reduce system complexity and minimize data transfer requirements, we can instead opt to use a single detector to perform homodyne detection (Fig. S1(b)). In this case, two copies of the signal $x_{1}$ and weight $A_{11}$ are sent into the detector, and a phase modulator applies a $\pi$-phase shift to the second copy of $x_{1}$, flipping its amplitude. The photocurrents $I_{+}$ and $I_{-}$ now appear separated in time rather than space, and are read out separately and may be subsequently subtracted. This technique can be applied over many channels and time steps, such that the full matrix-vector product $\vec{y}=A\vec{x}$ is computed (Fig. S1(c)). In this technique, half of the light is discarded and twice as many pulses are required, doubling the energy consumption and latency. However, since the same detector and beamsplitter are used for $I_{+}$ and $I_{-}$, one avoids the technical issues associated with beamsplitter imbalance and detector inhomogeneities.

In practice, optical fan-out of $\vec{x}$ and encoding of the weight-matrix $A$ by intensity and phase modulation can be performed on-chip. Optical phased arrays containing as many as 1024 grating antennas for outcoupling have been shown in a sub-cm2 CMOS chip, including 1192 phase and 168 amplitude modulators [84]. The signals can be interfered and imaged on the detector array using bulk optics. An example experimental setup is illustrated in Fig. S2.

The pulse trains for both $x_{i}$ and $A_{ij}$ must come from the same master laser, which can be either a pulsed laser or a continuous-wave laser that passes through an intensity modulator. For high-performance systems, moderate powers are required (e.g. assuming 0.1 pJ optical energy per detector, a 1-GHz $1000\times 1000$ optical GEMM will require $100\,{\rm mW}$ of laser power), so an optical amplifier may be needed. The path length difference must be less than both (1) the coherence length of the laser and (2) $c\tau_{\rm pulse}$. Although these conditions are not very stringent (assuming $\tau_{\rm pulse}\gtrsim 100\,{\rm ps}$), using a design with equal path lengths for data and weights is straightforward in practice.

S2 Aberration Management

Optical aberrations are a possible limit to scalability of our free-space homodyne optical neural network. In order to assess their impact, we used Zemax® ray-tracing software to model and optimize the lower branch of Fig. S2, which provides the weights $A_{ij}$. Specifically, on- and off-axis sources were propagated through an off-the-shelf achromatic lens (Thorlabs AC508-100-A) for collimation, a custom Cooke triplet of spherical lenses for focusing, and a cube beamsplitter. To achieve this design, illustrated in Fig. S3, we varied lens spacing and radii of curvature to minimize wavefront distortion due to Seidel aberrations across a $20\,{\rm mm}\times 20\,{\rm mm}$ field of view. To focus the beam in one dimension and collimate it in the other (required for the optical GEMM unit proposed in Fig. 6), the three circular lenses used for focusing can be replaced with cylindrical lenses with the same thicknesses and radii of curvature.

In the system described and modeled here, 1000 sources at $\lambda=635\,{\rm nm}$, each with a mode field diameter of $4\,{\rm\mu m}$ and separated by $13\,{\rm\mu m}$, are imaged onto a detector with one million $20\,{\rm\mu m}\times 20\,{\rm\mu m}$ pixels. The active area of each pixel is confined to the center $5\,{\rm\mu m}\times 5\,{\rm\mu m}$ to minimize crosstalk from neighboring pixels. As shown in Fig. S4, the Zemax® ray-traced aberrated spot sizes on the detector were found to be below the diffraction limit, and Zemax® “Physical Optics” simulations (based on scalar diffraction theory) confirm that we achieve small, near-Gaussian focused spots over the entire field of view. The substitution of circular lenses with cylindrical lenses does not lead to a significant decrease in lens performance, as the RMS spot size (in the dimension being focused) remains small.

The “Physical Optics” simulations were also used to estimate crosstalk, defined as undesired light intensity from adjacent pixels divided by signal intensity on the pixel’s active area. We found that >50% of the optical energy is detected in the active area, with $\sim$1% crosstalk from 8 nearest neighbors over the entire field of view. Interestingly, most of the crosstalk on the detector edges is due to a single neighbor, as the spot is aberrated asymmetrically.

Because of the phase coherence between the neighboring spots, constructive interference between a spot and its neighbors may make the actual crosstalk larger than the intensity figure quoted above. If all the light on the detector is assumed to be in phase (worst case), simulations show the crosstalk increases to $\sim$2%. If this level of crosstalk is too high for an application, it can be precalculated and compensated (to first order) in logic in the transmitter array. Based on these results, optical networks with up to $10^{6}$ channels should be achievable with reasonable signal fidelity. Further engineering improvements, e.g. apodizing the Airy disk with a Hann window, may lead to reduced crosstalk and larger system size limits.

S3 Shot Noise in Neural-Network Layers

S3.1 Matrix-Vector Multiply

In the time-encoded neural network, at each layer, a stream of data $\bar{x}_{i}$ is broadcast to the neurons of the subsequent layer. Each neuron is a homodyne detector that interferes this broadcast signal against the weight signal $\bar{A}_{ij}$. Assume perfect spatial and temporal mode-matching between input and weight signals, and normalize the quantities so that $|\bar{x}_{i}|^{2}$, $|\bar{A}_{ij}|^{2}$ correspond to the number of photons per pulse. If a pulse with amplitude $\bar{u}$ enters the detector, the output current is Poisson distributed:

[TABLE]

This has a mean of $|u|^{2}$ and a standard deviation of $|u|$. The homodyne signal at neuron $i$ is obtained by interfering the signals $\bar{x}_{j}$ and $\bar{A}_{ij}$ on a 50/50 beamsplitter and taking the difference between the photocurrents $y_{i}=\tfrac{1}{e}(Q_{i}^{(+)}-Q_{i}^{(-)})$. Each photocurrent $Q_{i}^{(\pm)}$ is the sum of many Poisson random variables, which is itself a Poisson process. In the useful limit of many photons per neuron (though not necessarily per MAC), this is approximately Gaussian:

[TABLE]

where $w_{i}^{(\pm)}\sim N(0,1)$ is a (normal) random variable. This term is responsible for the detector shot noise. The homodyne detector measures the difference between the photocurrents. The means subtract while the noises add in quadrature:

[TABLE]

In the final line of (S3), the sum over $j$ is implicit, and we have used the $L_{2}$ vector norms for $\bar{x}$ and the row-vectors of $\bar{A}_{i}$ to simplify the notation as follows: $\lVert{u}\rVert^{2}=\sum_{j}|u_{j}|^{2}$. As before, $w_{i}^{(\pm)}$ and $w_{i}$ are normally distributed as $N(0,1)$.

This output is sent through a nonlinear function (here a ReLU) to obtain the inputs to the next layer:

[TABLE]

Here $\bar{x}$, $\bar{x}^{\prime}$ and $\bar{A}$ are the physical inputs/outputs and weights, normalized to the single-photon level (so that $|\bar{x}_{i}|^{2}$ is the photon number in pulse $x_{i}$, etc.). These are related to the logical variables by scaling constants $\xi_{x}$, $\xi_{A}$:

[TABLE]

Because the ReLU is scale-free ($\mbox{ReLU}(cx)=c\,\mbox{ReLU}(x)$ for $c>0$), the relation between feature vectors at subsequent layers is:

[TABLE]

Set $\alpha=1/(2\xi_{A})$ to make the classical term match the desired relation $x^{\prime}_{i}=\mbox{ReLU}(A_{ij}x_{j})$. The quantum noise enters inversely in $\xi_{x},\xi_{A}$, which is expected because these quantities are related to the photon number, through the scaling of $(\bar{x},\bar{x}^{\prime},\bar{A})$.

The total number of photons involved in this layer is sum of the input photons and weight photons $n_{\rm tot}=n_{\rm tot}^{(x)}+n_{\rm tot}^{(A)}$, each given by:

[TABLE]

A more computationally relevant figure is the number of photons per MAC $n_{\rm mac}=n_{\rm tot}/NN^{\prime}$, which can also be divided into input and weight photons $n_{\rm mac}=n_{\rm mac}^{(x)}+n_{\rm mac}^{(A)}$:

[TABLE]

It is often appropriate to assume that the row vectors $A_{i}$ have approximately the same norm: $\lVert{A_{i}}\rVert\approx\lVert{A_{j}}\rVert$. In this case, we can replace $\lVert{A_{i}}\rVert\rightarrow\lVert{A}\rVert/\sqrt{N^{\prime}}$. Making this approximation and substituting $\xi_{x}$, $\xi_{A}$ from Eq. (S8), we find:

[TABLE]

For a fixed energy per MAC $n_{\rm mac}=n_{\rm mac}^{(x)}+n_{\rm mac}^{(A)}$, the SNR is maximized when $n_{\rm mac}^{(x)}=n_{\rm mac}^{(A)}=\tfrac{1}{2}n_{\rm mac}$, which gives the form used in the main text:

[TABLE]

S3.2 Matrix-Matrix Multiply

Fig. 6(a) of the main text introduced a modified scheme that multiplies two matrices rather than a matrix and a vector. This is advantageous when running the network on batches of data. It is also well-adapted to computing convolutions through the patching method (see main text). To compute the vector-vector product $C_{m\times n}=A_{m\times k}B_{k\times n}$, pixel $(i,j)$ receives the homodyne product of two pulse trains, encoding $A_{i,:}$ (the $i^{\rm th}$ row of $A$) and $B_{:,j}$ (the $j^{\rm th}$ column of $B$). The result is $C_{ij}=\sum_{l}A_{il}B_{lj}$, assuming all matrices are real.

As in Sec. S3.1, we normalize all quantities $\bar{A}_{ij}$, $\bar{B}_{ij}$, $\bar{C}_{ij}$ to photon number. Following Eq. (S3), the measured value of $\bar{C}_{ij}$ is photocurrent at detector $(i,j)$:

[TABLE]

The physical quantities are related to their logical values by a scaling factor: $\bar{A}=\xi_{A}A,\bar{B}=\xi_{B}B,\bar{C}=\xi_{C}C$. Scaling to logical quantities and setting $\xi_{C}=2\xi_{A}\xi_{B}$ (to scale $C$ to satisfy $C=AB$ in the classical limit) we obtain:

[TABLE]

The optical energy per MAC is $n_{\rm mac}=n_{\rm mac}^{(A)}+n_{\rm mac}^{(B)}$, where $n_{\rm mac}^{(A)}=\lVert{A}\rVert^{2}\xi_{A}^{2}/mk$, $n_{\rm mac}^{(B)}=\lVert{B}\rVert^{2}\xi_{B}^{2}/nk$. This allows us to trade $\xi_{A},\xi_{B}$ for $n_{\rm mac}^{(A)},n_{\rm mac}^{(B)}$ in Eq. (S12). Also, assuming that the rows of $A$ (and columns of $B$) have roughly the same norm, we can replace $\lVert{A_{i,:}}\rVert\rightarrow\lVert{A}\rVert/\sqrt{m}$ and $\lVert{B_{:,j}}\rVert\rightarrow\lVert{B}\rVert/\sqrt{n}$. Eq. (S12) then simplifies to:

[TABLE]

As before, subject to a fixed energy constraint on $n_{\rm mac}$, setting $n_{\rm mac}^{(A)}=n_{\rm mac}^{(B)}$ maximizes the SNR. Under this optimized choice of parameters, the input-output relation becomes:

[TABLE]

S4 Johnson Noise

Johnson-Nyquist noise will lead to thermal fluctuations in measured photocurrent, which imposes its own limit on energy consumption. However, unlike shot noise, this is hardware-dependent and can be mitigated by judicious detector design. This section modifies the analysis of Sec. S3.1 to include Johnson noise and determines the conditions under which it can be neglected.

Unlike shot noise, the amplitude of Johnson noise does not depend on the photocurrent. Thus we modify Eq. (S3) to include a constant noise term:

[TABLE]

Here $\langle\Delta n_{e}^{2}\rangle$ is the electron count variance due to thermal noise. Noise from the two detectors adds in quadrature, leading to the factor of 2. Making the substitutions (S4-S5, S8), and assuming $\lVert{A_{i}}\rVert\approx\lVert{A}\rVert/\sqrt{N^{\prime}}$ as above, we arrive at the formula:

[TABLE]

As before, we maximize the SNR subject to the constraint $n_{\rm mac}^{(A)}+n_{\rm mac}^{(x)}=n_{\rm mac}$. The optimal choice is the same found previously: $n_{\rm mac}^{(A)}=n_{\rm mac}^{(x)}=\tfrac{1}{2}n_{\rm mac}$. This leads to

[TABLE]

This differs from Eq. (S10) by the factor $\sqrt{1+2\langle\Delta n_{e}^{2}\rangle/Nn_{\rm mac}}$ in the noise term. Johnson noise is dominant when $Nn_{\rm mac}\ll 2\langle\Delta n_{e}^{2}\rangle$, and shot noise is dominant when $Nn_{\rm mac}\gg 2\langle\Delta n_{e}^{2}\rangle$. The value of $\langle\Delta n_{e}^{2}\rangle$ depends on the detector design. For typical optical receivers with amplifiers, it can be quite large. For femtofarad “receiverless” detectors being studied for on-chip interconnects, one can estimate $\langle\Delta n_{e}^{2}\rangle$ from the capacitor kTC formula $\langle\Delta n_{e}^{2}\rangle=kTC/e^{2}$ [70]. The ratio of Johnson noise to shot noise at the SQL is thus:

[TABLE]

where we have defined $C_{0}\equiv Nn_{\rm mac}e^{2}/2kT$. $C_{0}$ is a measure of how small a detector’s capacitance must be to achieve a quantum-limited (rather than thermally limited) performance, and depends on the neural network layer. Values for the networks studied in Sec. II are calculated in Table S1. From these figures, we see that $C_{0}$ is in the few-femtofarad range, sub-femtofarad receiverless detectors should be able to reach the SQL, while larger detectors will be limited by Johnson noise.

By going to cryogenic temperatures, it is possible to reduce the Johnson noise and thus the on-chip operating power. However, the noise-limited power scales linearly as $P\propto T_{c}$ ($T_{c}$ is the chip temperature), so the heat dissipated goes as $dQ_{c}\propto T_{c}$ and the entropy created in the chip $dS=dQ/T_{c}$ is temperature-independent. Since this entropy must ultimately be removed from the system, the heat dissipated by the cryostat (at room temperature) will be larger by at least $T_{h}/T_{c}$: $dQ_{h}\geq(T_{h}/T_{c})dQ_{c}$, where $T_{h}$ is room temperature (in practice the dissipation is much larger due to the low efficiency of most cryostats relative to a Carnot engine [88]). Since this lower limit to $dQ_{h}$ is independent of $T_{c}$, reducing Johnson noise by lowering the chip temperature does not confer any benefit in total energy efficiency.

S5 Landauer Limit Dependence on Architecture and Bit Precision

The Landauer Limit [71] sets a thermodynamic lower bound of $kT\log(2)$ for the energy of an irreversible logic gate with a single-bit output. For a more complex operation such as multiplication, the limit will be $(kT\log(2))\times G$, where $G$ is the gate count. The estimate of $G=10^{3}$ for a MAC, used in the main text, is roughly accurate for 32-bit precision. However, the trend in deep learning has been toward lower-precision arithmetic, which is faster and more energy-efficient [17]. Table S2 lists the gate counts and transistor counts of commonly used integer multpliers [72, 73, 89, 90, 91] at all common bit precisions.

Gate count is an algorithmic construct independent of the physical implementation. Any digital computer based on irreversible gates will satisfy the bound $E_{\rm op}\geq(kT\log(2))\times G$, regardless of how (or whether) the gates are implemented with CMOS transistors. Fig. S5 shows the Landauer limit for the multipliers listed in Table S2. Since multipliers require more gates than adders, the multiplier accounts for the majority of the gates for a MAC, and Fig. S5 is a rough lower bound for the Landauer limit to $E_{\rm mac}$.

Note that for all multiplier types the gate and transistor count scale quadratically with number of bits. By moving from 32-bit arithmetic to 8-bit arithmetic, the Landauer limit decreases by a factor of 16. For the most efficient multiplier under this metric (Wallace/Booth), it reaches just under 100 zJ/MAC ($10^{-19}$ J/MAC). This is significantly lower than the 3 aJ value quoted in the main text, but still slightly above the SQL obtained for the larger networks in Fig. 5 (main text), suggesting that even when comparing against low bit-precision digital architectures, the optical network can beat the Landauer limit.

In practice, digital designers use different multipliers depending on if they are optimizing for speed, area on chip, or energy consumption [91]. Wallace Tree based multipliers with Booth encoding have proved to be very energy efficient in modern designs and are common. Braun Multipliers are sometimes used for unsigned multiplication [90].