Open Access
Add to BibSonomy
Abstract
The Ising problem, a vital combinatorial optimization problem in various fields, is hard to solve by traditional Von Neumann computing architecture on a large scale. Thus, lots of application-specific physical architectures are reported, including quantum-based, electronics-based, and optical-based platforms. A Hopfield neural network combined with a simulated annealing algorithm is considered one of the effective approaches but is still limited by large resource consumption. Here, we propose to accelerate the Hopfield network on a photonic integrated circuit composed of the arrays of Mach–Zehnder interferometer. Our proposed Photonic Hopfield Neural Network (PHNN), utilizing the massively parallel operations and integrated circuit with ultrafast iteration rate, converges to a stable ground state solution with high probability. The average success probabilities for the MaxCut problem with a problem size of 100 and the Spin-glass problem with a problem size of 60 can both reach more than 80%. Moreover, our proposed architecture is inherently robust to the noise induced by the imperfect characteristics of components on chip.
© 2023 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement
1. Introduction
For the past few years, solving Ising problem [1], a simple and classical physical model, has received more and more attentions. Ising model was originally proposed to explain ferromagnetic phase transitions, and it has gradually been found to play an important role in many other fields. Furthermore, a number of classical combinatorial optimization problems can be polynomially reduced to an Ising problem [2], thus solving Ising problem becomes very significant [3]. However, like most of combinatorial optimization problems, the Ising problem is a Non-deterministic Polynomial-time (NP) problem, that is, the solution of Ising problem can only be verified in polynomial time. This means that large scale Ising problem is an exceedingly resource-consuming problem for traditional Von Neumann computing architecture and even can’t be solved along with the increase of the scale of problem [4].
Optical platforms have the advantages of ultra-high bandwidth, parallel signal processing with multi-degrees of freedom and low power consumption, making optical computing one of the candidates in “Post Moore Era” [5–9]. In this context, various architectures based on optical platforms have been proposed to solve Ising problem, including spatial-photonic Ising machines [10,11], coherent Ising machines [12,13], and so on. Due to the ability of processing super-scale problems, spatial-photonic Ising machines have been extensively studied. However, the coupling coefficient of this Ising machines, which is mapped with the Ising problem, is not completely arbitrary, limiting it to solve only specific problems. By employing FPGA, which can map arbitrary coupling coefficient into the measurement feedback, coherent Ising machines have overcome this problem. However, their success probabilities are not ideal, especially in experiment where many settings of the setup can't reach the ideal ones at the same time. For example, since the amplitude heterogeneity is not actively suppressed [14], only 21% success probability is obtained when solving Möbius Ladder problem with graph size 100 in Ref. [13]. Besides, comparing with the above scheme consisting of multiple bulky components, Ising machines fabricated on integrated platform are more compact and more stable, attracting lots of attentions. In Ref. [15,16], the key function, matrix-vector multiplication, is realized in optical circuit which greatly reduces the power consumption of algorithm. However, the ground state energy is found by utilizing the large noise approximation, leading to the solution not accurate enough especially for the coupling coefficient of the Ising problem with high precision. On the other hand, the output of network eventually converges to a stable distribution rather than a state, which mean that further statistic of the sampling results is needed to obtain the minimum solution, making the Ising machine not efficient enough.
Hopfield neural network (HNN) [17], a single-layer full feedback neural network, is one of cornerstone of modern recurrent neural network. Discrete HNN is initially designed with associative memory functions, which can be used for pattern recognition, error correction and so on. The associative memory function refers to the HNN for arbitrary input states converging to a stable result called an attractor, corresponding to the local or global optimum of Ising problems. Thanks to the capacity of discrete HNN to dynamically find the local or global minimum for the arbitrary Ising problems, the combination of discreate HNN with annealing system, which employs simulated annealing (SA) algorithm [18], can be used to solve Ising problems. This combination has been demonstrated on memristor arrays, where minimum solution can be found with high success probability and no mathematically approximation [19]. Here, to more efficiently and more accurately solve the Ising problem with high speed, we proposed a photonic Ising machine by employing Hopfield neural network with simulated annealing algorithm. Through a large number of simulations, we proved the feasibility of the structure. The average success probabilities of both MaxCut problems with N = 100 and spin glass problems with N = 60 are exceed 80% by setting appropriate annealing rates. Besides, we further show that it is robust to a certain amount of noise of optical devices thanks to the annealing process. For MaxCut problems with N = 20, the success probability has no evident change when σps = 0.03 and σdc = 0.015, which results in a worst-case phase shift of about 0.1 rad for phase shifter and 10% ratio shift for directional coupler, and the above error parameters can be easily achieved by current commercial factories.
2. Architecture and operating principle
2.1 Theory of HNN-SA
To solve the Ising problem, the HNN augmented with noise which imitate simulated annealing process (HNN-SA) is employed in our proposed PHNN. The topology of HNN-SA is shown in Fig. 1(a). There are N neurons, and the output of each neuron j (j = 1, 2, …, N) on the next time are decided by the Gaussian noise, excitation function and the outputs of all the neurons except itself. This network can be mathematically described as [20]:
Fig. 1. (a) The topology of HNN working in an asynchronous mode. (b) Energy descent processes for an Ising problem for three cases: direct iterative solution of HNN, global search of SA algorithm, and HNN introduced into annealing process. Here, the histogram of inset show the energy result of 100 repeated runs of HNN and HNN-SA, respectively, the two-dimensional block in the inset shows the weight matrix K formed by kij, and initial standard deviation, final standard deviation, and annealing rate p are set as 4, 0.2 and 0.999, respectively. (c) The schematic diagram that the solution jumps out of local optimum with the help of noise.
In the equations above, xi and xj are the state of ith and jth neuron, respectively, nj(t) is the additional noise of the output of neuron, bj is the threshold of jth neuron, which are shown in Fig. 1(a), sgn is the excitation function of network and kij is the weight between neuron i and j. The Eq. (3) shows the asynchronous working mode of HNN-SA, that is, HNN-SA runs with only one neuron state excited according to the above equations at each time step. This asynchronous working mode can be implemented by switch arrays as shown in Fig. 1(a). At the end of each cycle, to imitate the simulated annealing process, the Gaussian noise nj(t) is set to gradually decrease, which is realized in the theoretical model of discrete HNN-SA by decreasing the standard deviation φ(t) with an annealing rate p, that is φ(t + 1) = pφ(t).
The energy function of a HNN-SA, which is equal to the discrete HNN, can be defined as [20] according to Eqs. (1-3):
If the additional noise of the output of neuron nj(t) = 0, this will lead to a traditional discrete HNN. It has been proved for the discrete HNN that △E(t)=E(t + 1)-E(t) ≤ 0 [20]. To further verify the validity of this characteristic, we have solved an Ising problem with the number of spins N = 40 and nj(t) = 0 by using Eqs. (1)–(4), and have plotted the energy descent process as shown by the blue line in Fig. 1(b). Here, to simplify the process of proving, the weights kij are randomly set as 1 or 0, and thresholds of all the neurons are set as 0 (bj = 0, j = 1, 2, …, N). It can be seen that the energy curve decreases monotonically with the increasing number of cycles and finally converges to a fixed point. This characteristic of the discrete HNN makes it possible to solve Ising problem, but it always falls into local optimum. This can be found by the blue curve which has not decreased to the minimum (shown by the black dashed line). Thus, the noise in the HNN-SA is of great important which can introduce perturbations to the network and help the iteration to escape the local minimum as shown in Fig. 1(c). By gradually decreasing the noise with p = 0.999 to mimic the SA, the HNN-SA can achieve the global optimum as the red solid line shown in Fig. 1(b). Although the number of cycles for HNN-SA to get the global optimum (∼1800), is larger than that for the discrete HNN to get the local optimum (∼200). These lead to a conclusion that the HNN-SA, which can obtain the global optimum at the cost of searching time, can balance searching ability of the discrete HNN and the SA algorithm. Thus, we propose to introduce the HNN-SA as the theoretical model into our photonic Ising machine.
2.2 Architecture of PHNN
To construct the HNN-SA in photonic method, an architecture which can be integrated on a chip is shown in Fig. 2(a). It’s clear that the most energy-intensive task here is matrix-vector multiplication. Considering that the matrix-vector multiplication is carried out over the real-number field, our photonic computational architecture adopts a photonic integrated circuit which can recognize the symbol of the outputs [21]. This photonic integrated circuit is composed of MZI arrays, including input preparation layer, weight multiplication layer and coherent detection layer. Figure 2(b) is the partial diagram of photonic integrated circuit, and a MZI is composed of two directional couplers with the power splitting ratio of 0.5 and two phase shifters. Here, MZI is considered as an ideal model regardless of polarization state, which should be selected in the further when conducting actual experiment. The red MZIs form the input preparation layer, and the network composed of blue MZIs represents the weight multiplication layer. In the input preparation layer, the input laser is split to N + 1 ports with the optical signal of a same amplitude by adjusting the internal phase shifters, leading to the input with an all-one vector. This means the spins are all with the up states. Subsequently, we modulate the external phase shifters to obtain the vector xi(t). Here, except compensating the phase produced by the ahead splitters and internal phase shifters, the phase values of the phase shifters are modulated with only two types of phase values (0 rad and π rad) to represent the two spin states of the input vector xi(t) (+1 and -1). It should be noted that compensated phases need only one calculation, leading to no additional computing consumption during the iteration. It’s worth noting that the last output port in the input preparation layer is the reference light, which is used for coherent detection with the output in the following weight multiplication layer.
Fig. 2. (a) Operation principle of the PHNN. The spin states are encoded to the amplitude and phase of input laser by input preparation layer, passed through the on-chip photonic domain and the off-chip electron domain, and the result is recurrently fed back to the input preparation layer. (b) Partial topology of photon domain on chip, including input preparation layer and weight multiplication layer.
After the process by the input preparation layer, the optical signal xi(t) is fed into the following weight multiplication layer. Because the MZIs array is a passive network, the total output optical power of the network can’t be greater than the input optical power. Thus, before mapping weight matrix W to the matrix multiplication layer, we firstly normalize the weight matrix as W’=W/β, where β is the spectral norm of matrix W. After that, we decompose the matrix W’ as W’=UΣV† by using singular value decomposition, where U and V are unitary matrixes, Σ is a real-value diagonal matrix, and V† is the conjugate transpose of V. Here, V†, Σ and U are realized by three MZI arrays, respectively, using the method designed by Reck et al [22]. In this method, an N × N unitary matrix is decomposed into a product of a series of 2 × 2 unitary matrices, where an arbitrary 2 × 2 unitary matrix can be realized by a MZI on photonic integrated circuit.
Subsequently, only one output of the weight multiplication layer is randomly selected by using optical switch array which can be realized by using MZI units [23]. And then the selected output is transferred to electric signal by using coherent detection method assisted by the reference light. Here, the coherent detection can be realized by using an on-chip homodyning and balanced detectors [24]. Otherwise, the electrical signal is perturbed by a noise of Gaussian distribution with standard deviation φ(t) as shown in Fig. 2(a). In the further experiment, this noise can be added by using a Gaussian noise generator [25]. Then the perturbed signal is processed by an excitation function sgn() with threshold bj which can be realized by an electric hysteresis comparator. One element of the output is then substituted by the corresponding element xj(t) by storing the element in the register, and the resulted data xj(t + 1) is input into the input preparation layer by reading the register. By this method, we can realize feedback loop at high speed, and the computing speed is mainly limited by the modulated rate of input preparation layer, the detectors in coherent detection layer, and the off-chip electronic domain. It should be noted that although the decomposition algorithm is also a resource consumption outside of the computing task, it’s a once for all thing for arbitrary task of the Ising problems, which we usually require multiple runs to ensure the accuracy of solutions. Therefore, the additional resource consumption of decomposition algorithm becomes insignificant compared to repeated annealing process, especially with the increase of problem complexity.
3. Simulation analysis
In order to investigate the performance of our PHNN, we firstly use our PHNN to solve a Möbius Ladder problem, which is a relatively simple problem in Ising problem and belongs to the maximum cut (MaxCut) problem. It should be noted that for MaxCut problem, one need to find the cut of a given graph into two subsets with the largest number of their connecting weighted edges, represented by matrix W, and it is the same for Möbius Ladder problem but specially with circulant graphs. To simulate and analyze our proposed PHNN, firstly, we build the on-chip photonic domain as shown in Fig. 2(a) using the corresponding photonic devices simulated by Lumerical Interconnect software, and implement the off-chip electronic domain directly using the corresponding function. It should be noted that the interactions of the Möbius Ladder problem belong to {0, 1}, leading to the weight matrix of HNN-SA, K∈{0,-1} [15]. Here, the number of vertices in Möbius Ladder problem is set as N = 8 as the left figure shown in Fig. 3(a), the initial and terminal standard deviation [φ(t0) and φ(tmax)] of the additional Gaussian noise are set as 1 and 0.1, respectively, and the annealing rate is set as 0.95, leading to the maximum cycle approximating 45. To clearly demonstrate the evolution of spin states, we set all the 8 vertices with the spin up state (+1) as the left figure shown in Fig. 3(a). It can be found in the middle figure of Fig.3(a) that the spin states tend to be stable with the increase of cycles, and that the spin states are stable after iterating 14 times, leading to a more stable solution than Ref. [15]. The eventually stable spin states, shown in the right figure of Fig. 3(a), lead to a minimum energy E = -8 and cut value zMC =10 which is equal to the result calculated by using the “Biq Mac” solver [26]. Here, the cut value zMC is defined as [26]:
Fig. 3. The simulation result of the Möbius Ladder graph with N = 8 vertices. (a) An iterative process to find a ground state solution. The blue lattice represents spin up (x=+1) and the yellow lattice represents spin down (x = -1). (b, top) The result of the balanced photodetector, and the input of activation function netj in an iterative process. (b, bottom) The number of the selected channel corresponding to the process in (b, top). (c) Histograms of the obtained solutions in 100 runs for this graph.
Next, to further explore the performance of our proposed PHNN, MaxCut problems and spin-glass problems with randomly connected edges on a larger scale are considered. Since the runtime of the simulation by Lumerical Interconnect software is relatively long (Intel Xeon W-2133 CPU), we build a mathematical model of our system based on Python Neuroptica package [27] which can get the same results as using Lumerical Interconnect software. 100 times repeatedly solving processes are implemented based on Python for the instance in Fig. 3(a), and 97% success probability can be obtained. In a MaxCut problem with N = 60 and graph density of 50%, whose graph is shown in Fig. 4(a), our proposed PHNN can find the minimum energy E = -177 with the cycle at about 4.5 × 104 as the blue solid line shown in the upper panel of Fig. 4(b). When randomly generating 5 instances, including the instance shown in Fig. 4(a), the min(E) can all be well found and the success probability for 100 runs increasing with the increasing of the annealing rate is shown by the solid blue line in the below panel. Otherwise, since the iteration cycle is in inverse proportion to the logarithm of the annealing rate ($ \propto $1/log2p), the increasing of annealing rate leads to a largely increasing of iteration cycles as the solid pink line shown. Specially, when the annealing rate increases to 0.99993, the iteration cycles significantly increases to 4.5 × 104, leading to about 1.2 hour time consuming for simulation by the above referred electronic computer. Moreover, we investigate the relationship between success probability and problem size with the graph density at 50% as shown in Fig. 4(c). With the problem size increasing, we exponentially increase the iteration cycles to keep the average success probability above 80% for 100 runs. That is, when keeping a same average success probability, larger size of the problem is, more time is needed. In such architectures, the PHNN could thus find ground states of arbitrary Ising problems with graph size N∼100 and average success probability more than 80% within less than 18.5 microsecond (=1/FPHNN × cycles), while the time consuming for electronic computing is about 1.37 second (=1/F × Nite × cycles). Here, FPHNN is the clock frequency of proposed PHNN, which is set as 10 GHz. This requires the bandwidth of the integrated MZIs in the input preparation layer and that of the BPD both larger than 10GHz, and the analogue bandwidth of arbitrary waveform generator larger than 10GHz, too. F is the operating frequency of the electronic computing, which is selected as 2.7 GHz, Nite is the number of operations in one iteration, which is 20000, and cycles is about 18.4 × 104. In this condition, our PHNN can be ∼105 times faster than the electronic computing. Compared with other Ising machines, the clock frequency of our PHNN is 20 times faster than the Hopfield neural network based on memristor array (500MHz) [19]. In addition, our proposed architecture is similar to the integrated optical Ising machine based on PRIS algorithm [15], but our method has a more stable output and more accurate theoretical ground state solution.
Fig. 4. (a)A random graph of the MaxCut problem with a connection density of 50% and N = 60. (b) Average success probability and the corresponding annealing rate for 5 instances like instance (a) for varying number of cycles. Here, the φ(t0) and φ(tmax) are set as 5 and 0.2 respectively. (c) Average success probability and the number of cycles for 5 instances like (a) with problem size. Here, the φ(t0) is set as 1, 5, 5, 5, 8 for N = 4, 20, 40, 60, 100 respectively and φ(tmax) is set as 0.2. (d) A random graph of Ising Spin-glass problem with a density of 50% and N = 60. (e) Average success probability and the corresponding annealing rate for 5 instances like instance (d) for varying number of cycles. Here, the φ(t0) and φ(tmax) are set as 3 and 0.1 respectively. (f)Average success probability and the number of cycles for 5 instances like (d) with problem size. Here, the φ(t0) is set as 1, 1, 1, 2, 3 for N = 4, 10, 20, 40, 60 respectively and φ(tmax) is set as 0.1.
For the spin-glass problem which is more complex than the above MaxCut problem, our proposed PHNN can also fast find the min(E) with a large success probability. It should be noted that spin-glass model is a family of materials whose magnetic properties are particularly complex and its interactions Kij uniformly distributed in [−1, 1], which means that the Eq. (1) has more local minima. For an instance shown in Fig. 4(d), the min(E)= -130.4 can be find within 1.7 × 105 cycles of iteration as the upper panel and the solid yellow line in the below panel shown in Fig. 4(e). Moreover, the change trend of success probability with the annealing rate is the same as that of the MaxCut problem shown in Fig. 4(b), while the needed cycles, to obtain the same success probability, is about four times larger than that of the MaxCut problem. Similar performance can also be found in the condition when considering the influence of the size of spins as shown Fig. 4(f). Therefore, our proposed PHNN indeed can efficiently solve the Ising problem with higher speed than conventional electronic computer.
4. Imperfection analysis
In the above discussion, we assume the components are ideal, which is not consistent with the characteristics of the real chips fabricated by the current state of art. To investigate the degradation of the architecture due to imperfections of the components, we add Gaussian noise nps with the standard deviation σps into the shifted phase of phase shifter and ndc with the standard deviation σdc into the power splitting ratio of the directional couplers. In summary, the transmission matrix of a MZI with imperfection is expressed as:
Fig. 5. An example of phase shifters configuration with Max-Cut N = 20 in the case of (a) SVD with ordered singular values and (b) SVD with disordered singular values. With the increase of phase shifter error and directional coupler error, the decrease of average fidelity in two case: (c) SVD decomposition with ordered singular values and (d) disordered singular values. (e) For 5 MaxCut N = 20 instances, the relationship between success probability and error in the case of SVD decomposition with disordered singular values. (f) The relationship between average success probability of 5 MaxCut N = 20 instances for 100 runs and phase shifter error in the case of comprehensive imperfect hardware. Here, the φ(t0) and φ(tmax) are set as 0.2 and 0.01 respectively.
The matrix Mj and Msimj represent the jth column of target matrix and simulation matrix, and the operation ‘•’ means the product of two vectors. 5 randomly generated MaxCut problems of N = 20 are used in this section. Figure 5(c) and (d) describes the changing process of average fidelity of ordered and disordered singular values as the variances of two noises increase. The fidelity derived from ordered singular value tends to decline faster, and this result is similar to that reported previously [30]. Along the diagonal of Fig. 5(b), where SVD is with disordered singular values, we obtain the relationship between the success probability and component error (Fig.5(e)). The success probabilities have no evident change when σps = 0.03 rad and σdc = 0.015, which results in a worst-case phase shift of about 0.1 rad and 10% ratio shift of directional coupler. It should be noted that the current error of phase shift and directional coupler are both larger than the typical error (σps = 5 × 10−3rad and 2×ndc = 7%, respectively). As the noise continues to increase, the success probabilities decrease approximate linear instead of cliff. Furthermore, we investigate the decrease of average success probability under a more comprehensive imperfect model as shown in the Fig. 5(f). Here, the α (α∈[0.04,0.05]) and ndc (σdc = 0.015) of each MZI in the model represent the inherent imperfection of optical circuit including the inherent insertion loss and the splitting ratio error of couplers. As can be shown from the point that corresponding to σdc = 0.015 and σps = 0.03 in Fig. 5(e) and Fig. 5(f), respectively, the introduction of insert loss α result in reduced average success probability from 0.952 to 0.892. It’s worth noting that the imperfection of components can be further compensated by algorithms such as genetic algorithm, so that the fidelity of matrix and success probability in solving problem can be significantly improved.
5. Conclusion
We have demonstrated that HNN can be sufficiently accelerated in an integrated optical circuit composed of MZI arrays to find the ground states of arbitrary Ising problems. When solving two typical Ising problems, Max-cut problem with N = 100 and spin-glass problem with N = 60, our proposed PHNN can stably get the ground state with more than 80% average success probability. Moreover, the imperfection of the composed devices is also considered. After considering the typical error of phase shifter and directional coupler, which are induced by fabrication imperfection, digital to analog converter limited precision and the thermal effects, our PHNN still remain large success probabilities. It should be noted that once the mapping of the Ising problem is completed, the input spin can be modulated by binary adjustment of the external phase shifter in the input preparation layer. Therefore, the rate of optical matrix-vector multiplication, which is the heaviest resource consuming part of HNN, mainly depends on the modulation rate of phase shifter. By estimating, our proposed PHNN can be ∼105 times faster than the electronic computing when solving the problem with size of 100. Our proposed PHNN not only offers a new competitive and effective accelerator for the solving of Ising problem, but also shows the potential of MZI-based optical computing architecture in specific domain.
Funding
National Natural Science Foundation of China (62171055, 61821001, 62135009, 61705015, 61625104, 61971065); State Key Laboratory of Information Photonics and Optical Communications (IPOC2020ZT08, IPOC2020ZT03).
Disclosures
The authors declare no conflicts of interest.
Data availability
Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.
References
1. I. Ernst, “Beitrag zur Theorie des Ferromagnetismus,” Z. Phys. 31(1), 253–258 (1925). [CrossRef]
2. A. Lucas, “Ising formulations of many NP problems,” Front. Phys. 2(5), 1–15 (2014). [CrossRef]
3. F. Gloverand and G. Kochenberger, Handbook of Metaheuristics (Springer, 2006).
4. A. Wiegele, “Biq mac library—a collection of max-cut and quadratic 0-1 programming instances of medium size,” (2007), retrieved https://biqmac.aau.at/.
5. J. Liu, Q. Wu, X. Sui, Q. Chen, G. Gu, L. Wang, and S. Li, “Research progress in optical neural networks: theory, applications and developments,” Photonix 2(1), 5–39 (2021). [CrossRef]
6. C. Li, X. Zhang, J. Li, T. Fang, and X. Dong, “The challenges of modern computing and new opportunities for optics,” PhotoniX 2(1), 20 (2021). [CrossRef]
7. T. Zhang, J. Wang, Y. Dan, Y. Lanqiu, J. Dai, X. Han, X. Sun, and K. Xu, “Efficient training and design of photonic neural network through neuroevolution,” Opt. Express 27(26), 37150–37163 (2019). [CrossRef]
8. T. Zhang, J. Wang, Q. Liu, J. Zhou, J. Dai, X. Han, Y. Zhou, and K. Xu, “Efficient spectrum prediction and inverse design for plasmonic waveguide systems based on artificial neural networks,” Photonics Res. 7(3), 368–380 (2019). [CrossRef]
9. Y. Dan, Z. Fan, X. Sun, T. Zhang, and K. Xu, “All-type optical logic gates using plasmonic coding metamaterials and multi-objective optimization,” Opt. Express 30(7), 11633–11646 (2022). [CrossRef]
10. D. Pierangeli, G. Marcucci, and C. Conti, “Large-scale photonic Ising machine by spatial light modulation,” Phys. Rev. Lett. 122(21), 213902 (2019). [CrossRef]
11. D. Pierangeli, G. Marcucci, and C. Conti, “Adiabatic evolution on a spatial-photonic Ising machine,” Optica 7(11), 1535–1543 (2020). [CrossRef]
12. T. Inagaki, Y. Haribara, K. Igarashi, T. Sonobe, S. Tamate, T. Honjo, A. Marandi, P. L. McMahon, T. Umeki, K. Enbutsu, O. Tadanaga, H. Takenouchi, K. Aihara, K. I. Kawarabayashi, K. Inoue, S. Utsunomiya, and H. Takesue, “A coherent Ising machine for 2000-node optimization problems,” Science 354(6312), 603–606 (2016). [CrossRef]
13. P. L. McMahon, A. Marandi, Y. Haribara, R. Hamerly, C. Langrock, S. Tamate, T. Inagaki, H. Takesue, S. Utsunomiya, K. Aihara, R. L. Byer, M. M. Fejer, H. Mabuchi, and Y. Yamamoto, “A fully programmable 100-spin coherent Ising machine with all-to-all connections,” Science 354(6312), 614–617 (2016). [CrossRef]
14. Y. Yamamoto, K. Aihara, T. Leleu, K.-I. Kawarabayashi, S. Kako, M. Fejer, K. Inoue, and H. Takesue, “Coherent Ising machines—optical neural networks operating at the quantum limit,” npj Quantum Inf. 3(1), 49 (2017). [CrossRef]
15. C. Roques-Carmes, Y. Shen, C. Zanoci, M. Prabhu, F. Atieh, L. Jing, T. Dubček, C. Mao, M. R. Johnson, V. Čeperić, J. D. Joannopoulos, D. Englund, and M. Soljačić, “Heuristic recurrent algorithms for photonic Ising machines,” Nat. Commun. 11(1), 249 (2020). [CrossRef]
16. M. Prabhu, C. Roques-Carmes, Y. Shen, N. Harris, L. Jing, J. Carolan, R. Hamerly, T. Baehr-Jones, M. Hochberg, V. Čeperić, J. D. Joannopoulos, D. R. Englund, and M. Soljačić, “Accelerating recurrent Ising machines in photonic integrated circuits,” Optica 7(5), 551–558 (2020). [CrossRef]
17. J. J. Hopfield, “Neural networks and physical systems with emergent collective computational abilities,” Proc. Natl. Acad. Sci. U. S. A. 79(8), 2554–2558 (1982). [CrossRef]
18. S. Kirkpatrick, C. D. Gelatt Jr., and M. P. Vecchi, “Optimization by simulated annealing,” Science 220(4598), 671–680 (1983). [CrossRef]
19. F. X. Cai, S. Kumar, T. Van Vaerenbergh, X. Sheng, R. Liu, C. Li, Z. Liu, M. Foltin, S. M. Yu, Q. F. Xia, J. J. Yang, R. Beausoleil, W. D. Lu, and J. P. Strachan, “Power-efficient combinatorial optimization using intrinsic noise in memristor Hopfield neural networks,” Nat. Electron. 3(7), 409–418 (2020). [CrossRef]
20. R. Rojas, Neural Networks-A Systematic Introduction (Springer, 1996), pp. 337–371.
21. H. Zhang, M. Gu, X. D. Jiang, J. Thompson, H. Cai, S. Paesani, R. Santagati, A. Laing, Y. Zhang, M. H. Yung, Y. Z. Shi, F. K. Muhammad, G. Q. Lo, X. S. Luo, B. Dong, D. L. Kwong, L. C. Kwek, and A. Q. Liu, “An optical neural chip for implementing complex-valued neural network,” Nat. Commun. 12(1), 457 (2021). [CrossRef]
22. M. Reck, A. Zeilinger, H. J. Bernstein, and P. Bertani, “Experimental realization of any discrete unitary operator,” Phys. Rev. Lett. 73(1), 58–61 (1994). [CrossRef]
23. S. Nakamura, S. Takahashi, I. Ogura, J. Ushida, K. Kurata, T. Hino, H. Takeshita, A. Tajima, M.-B. Yu, and G.-Q. Lo, “High extinction ratio optical switching independently of temperature with silicon photonic 1 x 8 switch,” in OFC/NFOEC, (IEEE, 2012).
24. K. Kikuchi, “Fundamentals of coherent optical fiber communications,” J. Lightwave Technol. 34(1), 157–179 (2016). [CrossRef]
25. D. U. Lee, J. D. Villasenor, W. Luk, and P. H. W. Leong, “A hardware Gaussian noise generator using the Box-Muller method and its error analysis,” IEEE Trans. Comput. 55(6), 659–671 (2006). [CrossRef]
26. F. Rendl, G. Rinaldi, and A. Wiegele, “Solving Max-Cut to optimality by intersecting semidefinite and polyhedral relaxations,” Math. Program. 121(2), 307–335 (2010). [CrossRef]
27. B. Bartlett, “Neuroptica: an optical neural network simulator,” Github (2019), retrieved https://github.com/fancompute/neuroptica.
28. R. Shao, G. Zhang, and X. Gong, “Generalized robust training scheme using genetic algorithm for optical neural networks with imprecise components,” Photonics Res. 10(8), 1868–1876 (2022). [CrossRef]
29. H. Zhou, Y. Zhao, X. Wang, D. Gao, J. Dong, and X. Zhang, “Self-configuring and reconfigurable silicon photonic signal processor,” ACS Photonics 7(3), 792–799 (2020). [CrossRef]
30. M. Y. Fang, S. Manipatruni, C. Wierzynski, A. Khosrowshahi, and M. R. DeWeese, “Design of optical neural networks with component imprecisions,” Opt. Express 27(10), 14009–14029 (2019). [CrossRef]
