The processing and optimization of complex systems have rapidly advanced with computationally expensive tasks that require resource-intensive, high-performance computing hardware^[1]. This demand initiates intensive research on non-Von Neumann computing systems, such as the vector-matrix multiplier^[2], the photonic neural network^[3], the quantum computer^[4], and the photonic Ising machine^[5–13]. Among these Ising solvers, the spatial photonic Ising machine (SPIM), by encoding the spins as a phase matrix in spatial light modulators (SLMs), offers substantial benefits of high connectivity and speed in ground state search^[8,14–19]. This scheme is compatible with an Ising model with fully connected interactions or an equivalent quadratic unconstrained binary optimization problem due to its high connectivity and scalability^[20,21].

Noise in computing hardware exerts significant and non-negligible effects on system performance. In coherent Ising machines, noise acts as an error source, which may induce dynamics beyond the regime of Ising spins, consequently failing Hamiltonian evolution^[22]. Conversely, optimal levels of noise are shown to enhance the performance of SPIMs, particularly in addressing frustrating spin problems^[23]. Chip-level photonic Ising samplers benefit from noise that accelerates the ground state search and explores a wider phase space^[24]. Our previous work reported a similar result, where quadrature spatial photonic Ising machines with spontaneous noise surpassed the simulated results by over 30% when handling the large-scale Max-Cut problems^[25]. Additionally, noise has been a subject of early investigation in the realm of neural networks, where it has been established as a regularization term for the objective function contributing to the accelerated convergence of the network^[26,27]. Nevertheless, research on the mechanism of noise in Ising machines remains lacking, despite numerous proposed approaches for utilizing or mitigating noise.

Indeed, a smoothed analysis theory has been introduced by Spielman and Teng as a tool for understanding the impact of noise on analog calculation algorithms^[28]. Smoothed analysis theory offers a theoretical framework that combines the best of both worst-case and average-case analyses, shedding light on the behavior of algorithms in practical scenarios^[29]. Smoothed analysis has been effectively applied to different local search algorithms to solve a variety of classical non-deterministic polynomial (NP) problems, including the Max-Cut, the traveling salesman, and the binary tree problems^[29–31].

In this Letter, we aim to investigate the effect of noise on spatial photonic Ising machines with a large-scale spin based on the smoothed analysis theory. We demonstrate that strategically introducing digital noise can enhance the performance of SPIMs utilizing the fluid interface and particle inpainting (FLIP) algorithm for solving the Max-Cut problem. The experimental results indicate that the performance of the SPIM exhibits a 49% enhancement compared to the classical heuristic Sahni–Gonzales (SG) algorithm when subjected to blind noise levels. Moreover, this improvement can be further escalated to 73% after additional digital noise injection to the optimal noise level. Additionally, by combining the advantages of analog and digital computation, we propose an optoelectronic co-optimization method. The optoelectronic hybrid structure achieves comparable performance to the noise-injected SPIM at optimal noise levels while improving stability by thousands of times.

By considering the effects of noise in SPIMs or other analog hardware, conventional worst-case or average-case analyses fall short of capturing the actual performance of the algorithms. This is because noise can introduce perturbations that lead to significant deviations from the ideal solutions, rendering the analysis incomplete^[29]. Here, smoothed analysis theory provides a more comprehensive approach to understanding the behavior of algorithms in the presence of noise. The primary idea behind the smoothed analysis is so far based on the “potential function” approach, where the Hamiltonian function plays the role of the potential in SPIMs, $$H=\sum _{⟨l,k⟩}{\epsilon }_l{\epsilon }_k{x}_l{x}_k,$$(1)where the spin variable $x_l=\pm 1$ is encoded onto the wavefront phase of $0/\pi$, and the interaction coefficient $J_{l,k}$ is set by the amplitude modulation ${\epsilon }_l$ and ${\epsilon }_k$, as shown in Fig. 1(a).

If the interactions between spins are uniform ($J_{l,k}=1$) or strong among adjacent spins (ferromagnet), there exists a smoother energy landscape. However, in a given Ising problem, the actual state that corresponds to the distribution of Hamiltonian values tends to be complex and chaotic due to a weak correlation between the spatial scale and energy scale^[32]. Consequently, numerous conspicuous energy barriers arise that pose challenges for traversal. This constitutes the disordered terrain within the realm of physics, characterized by randomly distributed minima and a considerable count of energy barriers. Such circumstances present a formidable impediment in the pursuit of problem resolution [see Fig. 1(b)]. Noise can be viewed as a form of tepid adversarial perturbation, whereby it is neither entirely adversarial nor completely benign. The expected performance of the perturbation to an instance is defined as the smoothing complexity. Considering the intrinsic noise in SPIMs as Gaussian distributed, the variance measures the magnitude of smoothing complexity. Thus, the Hamiltonian in SPIMs has a smoothed form, $$\tilde{H}=\sum _{⟨l,k⟩}{\epsilon }_l{\epsilon }_k{x}_l{x}_k+\xi ,$$(2)where $\xi$ is a Gaussian random vector of variance ${\sigma }^2$. By systematically varying the parameter $\sigma$ over the range from zero to infinity, the smoothed analysis offers a valuable tool for interpolating between worst-case and average-case scenarios. For $\sigma =0$, the analysis reverts to the conventional worst-case scenario. As $\sigma$ increases, the influence of the random perturbation $\xi$ becomes predominant over the original $H$, leading to an average-case analysis. The theory enables the identification of smoothness levels beyond which an algorithm’s performance deteriorates significantly, helping to establish thresholds for noise tolerance in practical scenarios.

The numerical simulation is employed to illustrate how the energy landscape of a given problem changes with and without noise injection at different levels. To ensure that all spin states are traversed, simulations are conducted on a high-rank spin-glass model with 10 spins across varying levels of Gaussian noise. The interaction coefficients $J$ between arbitrary spins take values in the interval $[-1,1]$. This complex model can be efficiently solved through eigenvalue decomposition, which is subsequently performed on the space-photon Ising machine utilizing time-division multiplexing^[33]. We traverse 1024 kinds of spin states of the spin glass model through simulation and compute the system Hamiltonian with Eq. (3), which are subsequently subject to normalization, $$\tilde{H}=\sum _{⟨l,k⟩}{J}_{l,k}{x}_l{x}_k+\xi .$$(3)

Note that in order to simplify the calculation, we omit the negative sign in the expression of the standard Ising model. Therefore, a higher calculated value with Eq. (3) indicates a greater proximity to the ground state. In addition, since the interaction coefficient $J_{l,k}$ may be negative, it is still possible to compute negative values for certain spin states. The normalized Hamiltonian for all spin states is shown in Fig. 2(a). According to the distribution characteristics of Hamiltonian, we classify the spin states into four categories. The solution with Hamiltonian $<0$ is defined when the normalized Hamiltonian is less than zero. A normal solution is deemed if the Hamiltonian is less than 50% of the maximum value, which corresponds to the most probable spin states with a probability exceeding 50%. Spin states with a Hamiltonian ranging from 50% to 90% of the maximum are regarded as suboptimal solutions, and the optimal solution is defined as the spin state with a Hamiltonian that exceeds 90%. This tolerance of 10% is maintained because the value of the Hamiltonian will deviate somewhat after the noise injection.

Then, we inject different levels of Gaussian noise into the spin-glass model, with the variance of the noise measuring the smoothness level. Figure 2(b) displays the percentages of different kinds of solutions at four smoothness levels of 0, 0.001, 0.01, and 0.1, and Fig. 2(c) shows the energy landscapes at four smoothness levels of 0, 0.001, 0.01, and 0.1.

A strategy organizes the 1024 spin states of the 10-spin Ising model by expressing the ten spin values as binary codes and arranging them by their corresponding decimal magnitudes. Figure 2(c) presents a more intuitive representation of the energy landscape, with the spin states grouped into a $32\times 32$ matrix as a three-dimensional bar chart. Notably, this arrangement does not affect the Hamiltonian distribution. It is observed that over 90% of the spin states are normal or with Hamiltonian $<0$, while only approximately 2% are considered optimal. As the smoothness level increases, the suboptimal solutions start to be smoothed out of the original characteristics, gradually transitioning toward normal solutions. Consequently, the likelihood of becoming trapped in local optima diminishes. However, the features of optimal solutions can also be erased by the noise with excessively high smoothness levels. In such cases, the Ising system cannot stably stay at the optimal solution, and even the searched optimal solution no longer corresponds to the original Ising model. Based on the simulation analysis, a satisfactory smoothness level is targeted around the order of ${10}^{-2}$.

Then, we perform 100 ground state searches for each smoothness level, and the resulting probabilities are indicated in Fig. 2(d). The number of iterations is fixed at 100, which is significantly lower than the number required to traverse all spin states. In the absence of noise (smoothness level = 0), the system guarantees evolution towards the ground state, although there remains a probability greater than 50% of converging to a suboptimal solution. As the smoothness level increases, the probability of the system converging to the optimal solution gradually rises, ultimately surpassing that of the suboptimal solutions. The peak is observed at a smoothness level of 0.04, and once this threshold is surpassed, the success probability experiences a rapid decline. We also calculate the probability of searching for solutions with Hamiltonian $<0$ and normal solutions. It is evident that a lower level of smoothness significantly reduces the acceptance of both normal and bad solutions. Conversely, excessive smoothing leads to the retention of certain normal solutions and even solutions with Hamiltonian $<0$ because the features of the optimal solution are erased by the noise at this instance.

As a demonstration, we solve the Max-Cut problems using our noise-injected SPIM with the FLIP algorithm. The detailed algorithm flow is shown in Sec. 1 of the Supplement 1. For a given instance of the Max-Cut problem, we define the number of steps of the FLIP algorithm as the maximum number of local improvements possible for any initial cut and pivot rule. This represents the longest path in the FLIP algorithm’s transition graph and requires exponential time before reaching a local maximum in the worst-case scenario^[33]. Nevertheless, we are interested in the smoothed number of steps in the FLIP algorithm. Based on the framework of smoothed analysis, with a tiny amount of random noise, the worst-case scenario of FLIP taking an exponential number of steps is unlikely to occur^[34] because real-world examples frequently contain noise from various sources, such as rounding mistakes, measurement errors, and numerical imprecision. The experimental setup and the noise injection method are demonstrated in Sec. 2 of the Supplement 1. The instances of the Max-Cut problems consist of graphs of 20,736 nodes with random integer weights chosen from the interval $[-10,10]$. To assess the quality of the solution results, we will utilize the classical heuristic algorithm SG algorithm^[5,35], which will be executed on a high-performance computing system (Intel i9-13900 K, 5.8 GHz), as an evaluation criterion. The detailed SG algorithm flow is shown in Sec. 3 of Supplement 1.

Figures 3(a) and 3(b) illustrate the effect of the injected noise on sparse and fully connected Ising models in experiments, respectively. The scenario where the level of injected digital noise is zero is regarded as the baseline SPIM for the blind noise state. Even at blind noise levels, the mean values are improved by approximately 30% in the case of sparse connections and 49% in the case of full connections, compared to the SG algorithm. Furthermore, the Max-Cut value exhibits a trend of initially increasing and subsequently decreasing as the level of injected noise rises, aligning with both theoretical expectations and simulation outcomes. Within the gray region lies the range of noise levels where the injected digital noise can provide gain, which can be considered the threshold for noise tolerance. The optimal noise levels are observed at 0.04 and 0.05 [solid points in Figs. 3(a) and 3(b)], respectively, resulting in an increase in gain under noise to 46% and 55% on average. Particularly, the evolution of the cut value under blind noise and optimal noise levels are shown in Figs. 3(c) and 3(d), respectively. The light blue (orange) region represents the interval of the result distribution from the ten experiments, with the mean indicated by the dark blue (orange) line.

Note that while injecting noise may lead to gains, it also results in an unstable distribution and, consequently, a notable increase in the standard deviation of the results. In fact, when the noise level exceeds 0.1, the performance of the photonic Ising machine system cannot be guaranteed. Figure 3(e) illustrates two possible scenarios that occur when the cut value evolves under excessive noise: 1) it stops converging quickly and fails to find a good solution that outperforms the SG algorithm, and 2) it exhibits an instability in convergence and regresses to an inferior solution even after achieving a superior one. Therefore, the gain from the noise comes at the cost of reduced stability. To ensure accuracy, one viable solution is to increase the number of experiments conducted.

Furthermore, we expand our experimental scope to address the weighted Max-Cut problems across a range of graph densities of $[0.5,1.0]$. The threshold for noise tolerance for different graph densities and the performance gain of the SPIM are summarized in Sec. 4 of the Supplement 1. Here, we define gain as the percentage of the Max-Cut value searched by the SPIM that exceeds the results provided by the SG algorithm. The gain of the noise-injected SPIM is observed to be maximized by 73%, reflecting a 24% improvement over the blind noise scenario.

The above experimental results indicate that the relative level of noise in relation to the system Hamiltonian significantly determines the degree of smoothing in the energy landscape. When the smoothing is below the optimal noise threshold, the noise primarily affects lower energy levels of the Hamiltonian, leaving larger energy barriers unchanged, which prevents any improvement in system performance. Conversely, when the noise level is sufficiently high to smooth out certain barriers without altering the energy landscape of the optimal solution—approaching the optimal noise level—the system can effectively escape local optima while remaining close to the optimal solution, resulting in performance enhancement. However, for more complex problems, the distinction between the suboptimal and optimal solutions is often minimal. Once the noise begins to sufficiently smooth the suboptimal solution, the optimal solution becomes susceptible to similar effects of the noise. This phenomenon contributes to the significant variance observed in the experimental results. Our findings demonstrate that the smoothing effect of injected noise on the energy landscape of the Ising model is both global and continuous. Consequently, developing a method that can precisely control both the noise level and its timing is essential for enhancing the overall efficiency and stability of the system.

In the above analysis, the photonic Ising machines have a higher likelihood of producing superior outcomes than the numerical simulations on CPUs due to their inherent smoothing mechanism. However, it appears to be out of its depth in the face of finer energy differences. To be compatible with the existing noise, a system must fulfill two prerequisites: effectively leveraging the smoothing effect of the photonic Ising machine and minimizing the variance in the results.

Digital computation is more adept at handling intricate details. Thus, we propose a hybrid structure that combines photonic and electric calculations to enhance the searches of energy landscapes with multiple ground states, as shown in Fig. 4(a). We analyze the results obtained from the photonic Ising machine calculations and put them in the electronic Ising model to re-search, aiming to compensate for inaccuracies in the optical calculations. The detailed algorithm flow is shown in Sec. 5 of the Supplement 1.

To validate the correctness and feasibility of this method, we compare the Ising model, the photonic Ising machine, the noise-injected photonic Ising machine, and the optoelectronic hybrid Ising machine in solving the fully connected Max-Cut problem with 20,736 nodes. The evolution of the cut value during the searching process of different Ising machines is illustrated in Fig. 4(b). The photonic Ising machine experiences a significant improvement in cut value due to its inherent smoothing effect, resulting in a superior solution when compared to the Ising model. When operating at the optimal noise level, the photonic Ising machine maximizes this advantage while concurrently exhibiting an unstable distribution of results. Most strikingly, comparable outcomes can be achieved by electrically processing the results generated by the photonic Ising machine.

Furthermore, we employ optoelectronic co-optimization to address the weighted Max-Cut problems across varied graph densities, comparing it with Ising models and hardware. Figure 4(c) presents the maximum results yielded by different schemes. The data shows the following: 1) Optoelectronic co-optimization, involving the optimization of both the Ising model and the photonic Ising machine, demonstrates performance enhancements of up to 57% and 16% when compared to the simple electric or photonic computational structure. 2) The optoelectronic co-optimization method also demonstrates superior performance compared to the noise-injected SPIM. 3) Optoelectronic co-optimization significantly improves the system stability, reducing variance by several thousand times compared to the noise-injected approach.

The noise in optical Ising machines and analog devices is essential for their potential use in large-scale computational tasks. Considering that intrinsic noise cannot be completely eliminated, the ideal noise level cannot be predicted and controlled. Here, we incorporate smoothing parameters employing digital noise injection to provide a more realistic and comprehensive assessment of noise behavior based on the theory of smoothed analysis. Our numerical simulations and experimental findings verify that modest noise enhances the performance of the photonic Ising machine by up to 24%. The optimal noise level is not a constant value and may be closely linked to parameters, such as the type, density, and scale of the Ising problems. However, it is noteworthy that the photonic Ising machine exhibits robustness against noise and provides dependable results within the noise range of 0.04 to 0.07. Understanding the smoothed impact of noise on the SPIM system enables the exploration of alternative methods to enhance system smoothness while mitigating the unanticipated influence of noise. Therefore, by analyzing how digital and analog computation performs across various search scenarios, we develop an optoelectronic collaboration method that is resilient and capable of delivering reliable results that closely approximate the optimal noise level.

Our research is valuable in exploring the theory of smoothed analysis as a tool for understanding the impact of noise on analog calculations and their built-in local search algorithms. It can inspire researchers to analyze the impact of noise on hardware performance, establish noise tolerance thresholds, and ultimately design structures and algorithms that are better equipped to handle real-world scenarios where noise is inherent.