1. Introduction

Combinatorial optimization problems exist across diverse fields in science, industry, and society [1]. However, these are often NP-hard problems, and the computational cost to find the optimal solution grows exponentially with the problem’s size. Various algorithms have been developed to find optimal solutions, such as simulated annealing (SA) [2] and evolutionary computation [3]. Subsequently, equations based on the Ising model were proposed to solve combinatorial optimization problems efficiently. Ising machines are dedicated hardware devices that use physical pseudo spins to implement the Ising model and minimize the energy in the model [4]. The Hamiltonian of the Ising model without an external magnetic field is represented by

Ising machines are developed using various physical systems, which can solve combinatorial optimization problems quickly. Because the ability to represent the problem in an Ising model depends on the flexibility of the spin-spin interaction $J$, using a suitable physical system is vital for the developing practical machines. An example of a phenomenon used in such systems is the quantum mechanical effect, associated with superconducting circuits [5] and trapped ions [6]. These methods use quantum annealing [7] based on quantum fluctuation. SA has also been emulated in semiconductor integrated circuits, including CMOS annealing machines [8] and digital annealers [9]. Among implementations using physical phenomena, photonics-based implementation approaches are considered one of the most promising for addressing large-scale problems, because of the capabilities of light in parallel and high-speed processing and the robustness of computing systems [10]. One example is a coherent Ising machine [11,12], where pseudo spins are implemented with optical pulses generated by degenerate optical parametric oscillators [13]. Another is an integrated nanophotonic recurrent Ising sampler [14], where pseudo spins are implemented with coherent optical amplitudes.

A promising technique for controlling light on a large scale is spatial light modulation, (often used for computing), which exploits the parallel propagation characteristics of light [15–17]. A promising use of this optical strategy is the spatial-photonic Ising machine (SPIM) [18], which represents spins through the modulation of light waves using a spatial light modulator (SLM). Spin-spin interactions are achieved by superimposing light waves through free-space propagation. Compared with other physical implementation systems, SPIMs provide a simpler configuration and high scalability in handling the number of spins, because they use the propagation parallelism of light based on Fourier optics [19]. Because of these features, SPIMs have attracted much attention, and many systems and methods to enhance the SPIM have been considered. For example, annealing methods [20,21], a spin encoding method [22], an interaction model using the transmission matrix of the scattering medium [23], and that using nonlinear optical effects [24] have been proposed.

However, the SPIM scheme has limited flexibility in expressing interaction $J$ to solve combinatorial optimization problems; its primitive version can only solve problems with a rank-1 interaction matrix. This is a limitation for practical applications because only a few types of problems can be mapped. Several methods have been developed to remove this limitation, including (1) the quadrature photonic spatial Ising machine that introduces orthogonal phase modulation and external magnetic fields [25], (2) the phase-encoding and intensity detection Ising annealer, which uses computer holograms to perform calculations for each spin [26], and (3) a new computing model using the Gauge transformations implemented by wavelength multiplexing [27]. These methods have advantages for the representation of problems but their scalability in handling spins is reduced. We address this issue by proposing, a multi-component model of SPIM that uses time-division multiplexing [28]. In the multi-component model, the rank of the interaction matrix to be implemented is improved by decomposing the Ising Hamiltonian into several components that can be handled in the primitive SPIM scheme. For general combinatorial optimization problems, the coefficient weights of components are important because they affect the quality of the obtained solution and convergence speed in executing the Ising machine. For controlling the coefficients, each component of the Ising Hamiltonian is obtained by time-division multiplexing. Then the individual components multiplied with each coefficient weight are summed on a computer. Consequently, the overall Ising Hamiltonian could not be obtained in a batch, and all components must be acquired separately.

This paper presents a space-division multiplexed SPIM (SDM-SPIM), as an option of the system configurations to optically calculate a sum of multi-component Hamiltonian at a time while maintaining the high flexibility of the interaction matrix. In the space-division multiplexing scheme, the beams encoding each component are independently controlled to adjust the individual optical intensities, and thus the weight coefficients can be physically multiplied simultaneously. Then, the sum of the Ising Hamiltonian for each component can be obtained, by superimposing these beams. Furthermore, SDM-SPIM makes it possible to physically tune an optical parameter, a set of weight coefficients relating to constraint conditions of problems, and a part of the optimization process can be replaced by a physical process that dynamically tunes an optical parameter. This study aims to validate the method and its capability using physical parameter tuning, which is achieved by implementing an SPIM with spatial-division multiplexing. We constructed a proof-of-concept system and applied it to knapsack problems—combinatorial optimization problems with a constraint term. Furthermore, we analyzed the impact of physical parameters on the search characteristics of our method and investigated methods to enhance the search performance in the SDM-SPIM.

2. Method

2.1 SPIM with space division multiplexing

In SPIMs, the Ising Hamiltonian is optically calculated by phase modulation of an amplitude-distributed light. The obtained Ising Hamiltonian is evaluated and used to update the next spin state. The spin $\sigma _j$ is encoded into phase modulation $\phi _j \in \{0,\pi \}$ by an SLM, where $\sigma _j = \exp (\imath \phi _j)=\pm 1$. When the amplitude distribution of light is $\xi _j$, the element $J_{j h}$ of the interaction matrix $J$ is proportional to $\xi _j$ and $\xi _h$. The Ising Hamiltonian in SPIMs can be expressed as follows [18].

The Ising model with such a Hamiltonian is known as the Mattis model [29]. The interaction matrix $J$ is limited to the form $J\propto \boldsymbol{\xi } \boldsymbol{\xi }^T$, where $\boldsymbol{\xi }=(\xi _1,\ldots,\xi _N)^T$, and SPIMs can accommodate only a real symmetric matrix of rank 1 as the interaction matrix.

This study adopts a multi-component model based on the linear combination of the Mattis model to represent arbitrary Hamiltonians [28]. The Ising Hamiltonian is represented by

We consider an implementation of the model of Eq. (3) with spatial multiplexing; the concept of an SDM-SPIM is depicted in Fig. 1. Amplitude distributions $\xi ^{(k)}_{j}$, which are mutually incoherent for different $k$, are synthesized by an amplitude modulator with spatial multiplexing. The phase-modulated light after the SLM is Fourier transformed and detected as the optical intensity distribution $\boldsymbol{I}^{(k)}(x)$, where $x$ is the position of pixels. $\boldsymbol{I}^{(k)}$ is the optical intensity after being physically tuned with a coefficient parameter $\alpha ^{(k)}$ using an optical intensity modulator such as a neutral density (ND) filter. The total optical intensity $\boldsymbol{I}^{\text {(total)}}$ of the spatial $K$-multiplexed light captured by the image sensor is obtained as a summation $\boldsymbol{I}^{\text {(total)}}=\boldsymbol{I}^{(1)} +\boldsymbol{I}^{(2)}+\cdots +\boldsymbol{I}^{(K)}$. If the coefficients composing the Ising Hamiltonian $\alpha ^{(k)}$ have the same sign, the result can be obtained in a single shot. Because the intensity is a non-negative value, if the coefficients contain both positive and negative signs, the intensity should be obtained for each sign. In this case, the sum of intensities $\boldsymbol{I}^{\text {(total)}}_+$ and $\boldsymbol{I}^{\text {(total)}}_-$ of the terms with positive and negative coefficients, respectively, are imaged separately and processed electronically. The Ising Hamiltonian can be obtained by capturing the images twice at the maximum. Therefore, the computation time in an SDM-SPIM is independent of the rank $K$ of the interaction matrix.

 

Fig. 1. Schematic of a spatial-photonic Ising machine by physically tunable space-division multiplexing (SDM-SPIM).

Download Full Size | PDF

In the primitive SPIM, the Ising Hamiltonian is expressed as $\mathcal {H} \propto \|\boldsymbol{I}-\boldsymbol{I_T}\|$, where the target intensity $\boldsymbol{I_T}$ is a delta function [18]. The total optical intensity is normalized by the total power. With space-division multiplexing, the normalization is represented as $\boldsymbol{I}^{\text {(total)}}_{\text {norm}}=\frac {\boldsymbol{I}^{(1)}(x)+ \boldsymbol{I}^{(2)}(x)+\cdots +\boldsymbol{I}^{(K)}(x)}{\sqrt {\int \left (\boldsymbol{I}^{(1)}(x)+ \boldsymbol{I}^{(2)}(x)+\cdots +\boldsymbol{I}^{(K)}(x)\right )^2\text {d}x}}$. Then, the Ising Hamiltonian is obtained in the SDM-SPIM as follows.

2.2 Knapsack problems

We used knapsack problems in experiments using the SDM-SPIM. The knapsack problems are combinatorial optimization problems to find the subset of items that maximizes the total value by satisfying each knapsack’s weight limit. This problem cannot be transformed to the Ising model with a rank 1 interaction matrix. The 0–1 knapsack problem with integer weights is formulated as follows:

Equation (10) actually contains an additional constant term, which does not essentially change the Ising Hamiltonian. The coefficients $A$ and $B$ must satisfy the following condition [1], to obtain a solution that adheres to the constraints:

At the SA process in the experiment, the Metropolis algorithm was used to update the spin states with an updating probability $P(\sigma )=\text {min}\left (1,\exp \left (-\frac {\Delta \mathcal {H}(\sigma )}{T}\right )\right )$ depending on the temperature parameter $T$, where $\Delta \mathcal {H}(\sigma )$ is the difference between the Ising Hamiltonian of the state $\sigma$ and the reference.

2.3 Optical setup

We constructed an SDM-SPIM system and applied it to knapsack problems. Figure 2 illustrates the optical setup of the two-space-division-multiplexed SPIM. In this system, two independent optical intensity distributions are required, and thus, two He-Ne lasers (wavelength: 632.8 nm) were used as light sources that are mutually incoherent. The individual beams pass through an expander and enter a spatial light modulator SLM1 (Holoeye, LC2012, pixel size $1024\times 768$, pixel pitch $36\ \mu \mathrm {m}$). SLM1 is segmented into two regions to provide an amplitude distribution $\boldsymbol{\xi ^{(k)}}$ for each beam. An ND filter is placed in the optical path of beam 2 to adjust the optical intensity ratio between the beams and control the contribution of the constraint and objective terms in the Hamiltonian calculation. Without the filter, the intensity ratio of the two terms, i.e., the coefficient ratio $\beta \left (=\frac {B}{A}\right )$, was 4. The beams modulated by SLM1 are coaxially superimposed by a mirror and beam splitter (BS1), and their phases are modulated by an SLM2 (Holoeye, Pluto-2, pixel size $1920\times 1080$, pixel pitch $8\ \mu \mathrm {m}$). The planes at SLM1 and SLM2 have a conjugate relationship given the $4f$ system with two lenses L2 ($f=150\ \mathrm {mm}$). Because the size of the corresponding modulation areas between SLM1 and SLM2 must be matched, the modulation area representing each spin (spin pixel size) was set to a common multiple of the SLMs’ pixel sizes. In this experiment, the spin pixel size was set to $360\times 360\ \mu \mathrm {m}^2$ ($10\times 10$ pixels for SLM1 and $45\times 45$ pixels for SLM2). The same phase modulation is applied to beams 1 and 2, and the beams are Fourier transformed by lens L3 ($f=300\ \mathrm {mm}$) and captured by a CMOS image sensor (PointGray Research, Grasshopper GS3-U3-32S4, pixel pitch $3.45\ \mu \mathrm {m}$, 16-bit monochrome grayscale). The optical intensity $I^{\text {(total)}}$ is obtained by cropping the image to $10\times 10$ pixels at the center. In the experiment, we used $\left |\xi ^{(k)}_{j}\right |$ as the amplitude with flipping the corresponding spin $\sigma _j$ in case that the amplitude modulation $\xi ^{(k)}_{j}$ is negative.

 

Fig. 2. Optical setup in the experiment.

Download Full Size | PDF

3. Result

3.1 Demonstration using knapsack problems

We investigated the system behavior using a knapsack problem with 13 items (Appendix A: Problem 1) to demonstrate the SDM-SPIM. In Problem 1, the weight constraint is $\mathcal {W}\leq 80$. The coefficient ratio was chosen to be $\beta =0.04$ (ND filter transmittance was 1%) per the requirements of Eq. (15). The required amplitude corresponding to the fixed spin was 45, but it was divided into three parts of 15 each considering the dynamic range of the modulation amount. Consequently, the total number of spins at execution was $N=17$. The initial temperature in SA was experimentally determined and the rate of decrease was set to a constant rate. The solution search was executed with 3000 iterations, and we solved the problem 50 times with randomly determined initial conditions.

Figures 3 (a) and 3 (b) illustrate histograms of the total weight and total value of the solutions obtained. The results demonstrate that the optimal solution (total value: 95) was obtained. In the demonstration experiment, the final solution was determined as the sample with the highest value under the constraint weight among all the explored samples. Figure 3 (c) illustrates the sample distributions obtained during iterations. The horizontal axis indicates the total weight of each sample, and the vertical axis indicates the total value of each sample. As the iteration progresses, the search area converges to an area around the optimal solution. Figure 3 (d) illustrates an example of the time course of the Ising Hamiltonian with iteration. As illustrated at the end of the iteration in Figs. 3 (c) and 3 (d), the near-ground state of Ising Hamiltonian obtained in the experiment certainly corresponds to samples of near-optimal solution.

 

Fig. 3. Result using the 13-item knapsack problem. (a) Histogram of the total weight for obtained solutions. (b) Histogram of the total value for obtained solutions. (c) Example of the distribution of explored samples. (d) Example of the time evolution of Ising Hamiltonian values.

Download Full Size | PDF

3.2 Optical parameter tuning

The coefficient ratio, $\beta =\frac {A}{B}$, of the constraint term to the objective term should be set appropriately, to obtain the optimal solution of a knapsack problem. In the SDM-SPIM, this ratio can be controlled by the intensity of the light corresponding to each term. We solved a four-item knapsack problem (Problem 2) to investigate the impact on the property of the solution search when changing the transmittance of the ND filter in the optical path of the laser 2. The number of the auxiliary variable for Problem 2 was set to $m=4$, and the total spin number was $N=8$ to check that the constraint is still not violated by adding a sizable auxiliary spin. The temperature parameter in SA was set to decrease at a similar constant rate for all ND filters.

Figure 4 (a) illustrates the number of explored samples that violated the weight constraint throughout the iterations, and Fig. 4 (b) is the histogram of the total weight for each solution obtained with the ratio $\beta =0.4$ and $0.01$ (ND filter transmittances were 10% and 0.25%). We solved each problem 50 times. The total iteration was 300. When the optical intensity representing the objective term is decreased, the search is more likely to be run while adhering to the constraint. In contrast, when the optical intensity is increased, the search is more likely to be run in violation of the constraint. In the result, when $\beta =0.01$ (ND: 0.25%), no solution violated the weight constraint. In contrast, when $\beta =0.4$ (ND: 10%), many solutions were obtained that violated the constraint. The experimental result illustrates that the distribution of the obtained solutions changed with the ND filter. This indicates that the solution-search space can be manipulated by tuning the optical parameter. An effective constraint condition should be $\beta \leq \frac {1}{13}=0.077$, because $\text {max} v_i=13$ in Problem 2, according to Eq. (15). The result in Fig. 4 is consistent with this.

 

Fig. 4. Result with ND filter tuning. (a) Frequency of samples that violate the weight constraint in the iterations. (b) Histogram of the total value for obtained solutions with 10% and 0.25% ND filters.

Download Full Size | PDF

In Fig. 4, we obtained results using a fixed ND filter transmission during iterations. However, if the ND filter could be changed for each iteration, it would be possible to change the search trend in steps. Based on this idea, we investigated a new annealing method, dynamic coefficient search (DCS), where the coefficient ratio is changed step-by-step during iteration. Figure 5 (a) provides an overview of the DCS process. Varying the coefficient ratio in steps during the iteration changes the value of the Ising Hamiltonian for the same spin state. The spin state with the minimum Ising Hamiltonian obtained in the search using the same coefficient ratio is adopted as the initial state for the subsequent coefficient ratio. The spin state is updated by annealing at a constant temperature. By repeating this process, the optimal solution is obtained. Changes in coefficient ratios can be achieved in the optical system by, for example, mechanically changing the ND filters or adjusting a laser with variable emission intensity.

 

Fig. 5. Result of numerical simulation of the Dynamic Coefficient Search (DCS). (a) Schematic of the DCS. (b) Histogram of the total value for solutions obtained by the DCS and the SA. (c) Comparison of frequency of the optimal solution obtained between the DCS and the SA. (d) Example of time evolution of Ising Hamiltonian values obtained with the DCS and the SA.

Download Full Size | PDF

We compared the SA and the DCS by solving a knapsack problem (Problem 3) with a numerical simulation. The coefficient in SA was fixed at $\beta =0.05$. The coefficient in the DCS was varied in the order $\beta =2, 1, 0.8, 0.5, 0.1, 0.05$ at each step, and the annealing temperature was kept constant.

Figure 5 (b) illustrates the histogram of the total value for solutions obtained in 1000 executions. The total number of iterations was set to 600, with 100 iterations for each coefficient in the DCS. The results illustrate that DCS is more likely than SA to obtain solutions closer to the optimal solution. Figure 5 (c) illustrates the frequency of optimal solutions obtained out of 1000 runs when the total number of iterations is varied. The trend remained unchanged in this result, with the DCS demonstrating superior properties to the SA. A possible reason for the improved result is the difference in the convergence process of the solutions. For example, Fig. 5 (d) illustrates the time evolution of the Ising Hamiltonian for SA and the DCS. In the knapsack problem, the Ising Hamiltonian of the constraint term converges to 0 as the knapsack weight approaches its limit: the smaller the objective term, the higher the value of the selected items. Based on Eq. (15), a solution with a smaller Ising Hamiltonian may not necessarily be superior because it may violate the constraints. In the SA with the fixed coefficients, the constraint term is dominant from the initial iterations, and the search is preferentially focused on solutions that satisfy the constraints. In contrast, in the DCS, the constraint is weak in the initial iterations, enabling it to search for solutions with a high total value while violating the constraint. This property may help to find the optimal solution quickly.

4. Discussion

Knapsack problems are a specific type of low-rank problem that can be efficiently represented as a multi-component form with rank $K=3$, as in Eqs. (10)–(14). While the multi-component model is particularly efficient for low-rank Ising Hamiltonians, it is also capable of dealing with any Ising Hamiltonian [28]. The matrix of the Ising Hamiltonian is real-symmetric, so that it can be converted to Mattis-type matrices by spectral decomposition, which is mainly based on the calculation of eigenvalues and corresponding vectors. Various algorithms exist for this purpose, and it is known that the typical computational cost is $O(N^3)$, which grows only polynomially with $N$ [30]. In particular, this can be reduced for small-rank matrices using, for example, the Krylov subspace method [31]. As a result, the cost of matrix conversion is less than that of the annealing process when dealing with large problems.

In this experiment, the spin pixel size was the same for the two SLMs. However, to increase the number of spins that can be handled, smaller spin pixel sizes should be realized by matching the pixel size of SLM1 to that of SLM2. A simple method is the use of magnified imaging of the SLM1 plane to the SLM2 plane. By using an appropriate magnification, the system with different spin pixel sizes can be constructed.

The SDM-SPIM increases the flexibility of the interactions by extending the interaction matrix with spatial resources by acquiring the Ising Hamiltonians in a single or double shot. In the demonstration, we applied three amplitude distributions that are spatially and temporally multiplexed to represent a rank-3 interaction matrix. At most two time-division processes for individual signs are sufficient. Thus, the time required for time-division processing is constant, regardless of the rank of the interaction. In addition, it will be possible to obtain the Ising Hamiltonian by a single shot utilizing an appropriate spatial division method. For example, by adding a bias phase which works as a blazed grating to a phase-modulating SLM, the focusing position is shifted on the image sensor. When a phase-modulating SLM is divided into two regions for spins of positive and negative terms, and different bias phases are added, all components of the Ising Hamiltonian can be obtained in a single shot by taking two intensity distributions that are focused at different positions, respectively.

The multi-component model works efficiently, in particular for low-rank combinatorial optimization problems. Thus, for low-rank problems, the system implemented using spatial multiplexing can work well even when using a few optical components. In this study, the optical system was constructed using an SLM, mirrors, and BSs in the multiplexing part. However, because of the limitation of the physical size of the elements, the beam illumination and combining part should be made more compact to increase the multiplexing number. Dynamic control of an amplitude-modulating SLM is not needed and it can be replaced by a static device. By assigning such devices to individual optical paths, the optical setup is expressed to be constructed more freely with beam splitters or other general devices. Note that dividing the region of a phase-modulating SLM is unnecessary, so the number of its pixels is not a limiting factor of the multiplexing number.

The DCS considered in this paper seems advantageous not only because of its potential to improve search efficiency but also because it can extend the physically implementable part of the system. Part of the DCS process can be physically achieved by adjusting the ratio of the coefficient of each term. The SDM-SPIM setup enables the coefficients to be adjusted simply by adjusting the ratio of the intensity of the encoded light. In the current system, the most straightforward method of dynamically changing the light intensity is manually replacing the ND filter sequentially. However, this method is not the best method because it leads to changes that are small but affect results in the light patterns due to replacement accuracy. We consider that a system with a laser whose intensity is electronically tunable or an intensity modulator such as an acoustic optical modulator (AOM) would be appropriate for stable tuning the intensity. Simulated fluctuation in annealing can also be replaced using noise fluctuation during imaging [20]. Other parts of the system, such as the evaluation of the Ising Hamiltonian and control of the devices will still rely on the computer as before. With the SDM-SPIM setup using intensity-variable lasers and the DCS, the Ising Hamiltonian calculation and electrical portion of the optimization process is physically archivable.

We can also consider other types of problems for which applying the DCS is more effective. Although knapsack problems have only one constraint term, the DCS has the potential to select constraints preferentially to be kept in a problem with multiple constraints. This potential is expected to be valuable in obtaining approximate solutions sufficient for practical use.

5. Conclusion

We presented a SPIM system with physically tunable space-division multiplexing as one of the implementing methods to physically calculate the weighted sum of the Ising Hamiltonian components at a time in a multi-component computing model. We experimentally demonstrated the validity of the SDM-SPIM by solving knapsack problems with a constraint term using the system. Another result demonstrates that the solution distribution obtained by the SDM-SPIM can be changed by controlling the intensity ratio of the encoded light using an ND filter placed in the optical path. Based on this property, we considered a new DCS algorithm as a solution-search method and compared it with the SA by simulation. The results suggest the ability to improve the performance of the SDM-SPIM by dynamically tuning optical parameters. Our method computes the Ising Hamiltonian in a single or double shot, saving computation time for high-rank problems. The SDM-SPIM that can handle a large number of spins will enable the solving of various combinatorial optimization problems at high speed by increasing the number of multiplexes.