Photonic restricted Boltzmann machine for content generation tasks

TL;DR

This paper introduces a photonic restricted Boltzmann machine (PRBM) that uses photonic computing to significantly accelerate Gibbs sampling, enabling efficient and scalable content generation for AI applications.

Contribution

The paper presents a novel PRBM architecture that reduces Gibbs sampling complexity to O(1) and eliminates matrix storage, advancing photonic AI hardware capabilities.

Findings

Photonic Gibbs sampling closely matches theoretical phase transition temperatures.

PRBM effectively generates and restores diverse content including images and sequences.

The approach demonstrates scalability and lower training costs for large-scale models.

Abstract

The restricted Boltzmann machine (RBM) is a neural network based on the Ising model, well known for its ability to learn probability distributions and stochastically generate new content. However, the high computational cost of Gibbs sampling in content generation tasks imposes significant bottlenecks on electronic implementations. Here, we propose a photonic restricted Boltzmann machine (PRBM) that leverages photonic computing to accelerate Gibbs sampling, enabling efficient content generation. By introducing an efficient encoding method, the PRBM eliminates the need for computationally intensive matrix decomposition and reduces the computational complexity of Gibbs sampling from $O(N)$ to $O(1)$. Moreover, its non-Von Neumann photonic computing architecture circumvents the memory storage of interaction matrices, providing substantial advantages for large-scale RBMs. We experimentally validate the photonic-accelerated Gibbs sampling by simulating a two-dimensional Ising model, where the observed phase transition temperature closely matches the theoretical predictions. Beyond physics-inspired tasks, the PRBM demonstrates robust capabilities in generating and restoring diverse content, including images and temporal sequences, even in the presence of noise and aberrations. The scalability and reduced training cost of the PRBM framework underscore its potential as a promising pathway for advancing photonic computing in generative artificial intelligence.

Equations

$$\begin{array}[]{ll}\text{Region}\;{\rm I}:&\{\phi_{k}^{m,n}\}^{\rm I}={h_{k}}\frac{\pi}{2}\\ \text{Region}\;{\rm II}:&\{\phi_{i,k}^{m,n}\}^{{\rm II}}={v_{i}}\frac{\pi}{2}+{(-1)^{m+n}}{\alpha_{ik}}\\ \text{Region}\;{\rm III}:&\{\phi_{k}^{m,n}\}^{{\rm III}}=s\frac{\pi}{2}+{(-1)^{m+n}}{\beta_{k}}.\end{array}$$

$$\Delta H_k=(H_0-H_k)/2=-2\sum _{j=1}^{N_v}{v}_j{W}_{jk}-2b_k.$$

$$E_k(x,y)=A_0[G_{+}\cdot H+G_{-}\cdot (1-H)]\cdot \mathrm{Rec}{\mathrm{t}}_{N_{x}W}(x)$$

$$G_{\pm }$$

$$=i{h}_k\cdot \mathrm{Rec}{\mathrm{t}}_{N_{y}W}\left(y-y_1\right)+i\sum _{j=1}^{N_v}{e}^{\pm i{\alpha }_{jk}}{v}_j\cdot \mathrm{Rec}{\mathrm{t}}_{N_{y}W}\left(y-y_{j+1}\right)+i{e}^{\pm i{\beta }_k}s\cdot \mathrm{Rec}{\mathrm{t}}_{N_{y}W}\left(y-y_{N_v+2}\right)$$

$${{\mathop{\rm Rect}\nolimits}_{{N_{y}}W}}(y)={\rm{Rect}}(\frac{y}{N_{y}W})=\left\{{\begin{array}[]{*{20}{l}}{1},&{|y|\leq{N_{y}}W/2}\\ 0,&{|y|>{N_{y}}W/2}\end{array}}\right.$$

$$\begin{array}[]{c}H=\sum\limits_{m,n=-\infty}^{\infty}{[\delta(x-2mW,y-2nW)+\delta(x-(2m+1)W,y-(2n+1)W)]}\otimes{\rm{Rec}}{{\rm{t}}_{W}}(x,y)\end{array}$$

$${\tilde{E}}_k(k_x,k_y)=\frac{1}{i\lambda f}\frac{A_0}{4}(P_{+}\otimes Q_1+P_{-}\otimes Q_2)$$

$$P_{\pm }=\sum _{j=1}^{N_v+2}{e}^{\pm i{\varphi }_{j,k}}\cdot \sin {\mathrm{c}}_{N_{y}W}\left(k_y\right)e^{i{k}_y{y}_j}\cdot (N_{y}W)\cdot \sin {\mathrm{c}}_{N_{x}W}(k_x)\cdot (N_x\mathrm{W})$$

$$Q_1=\sum _{m,n=-\infty }^{\infty }\left(1+{(-1)}^{m+n}\right)\delta \left(k_x-m\frac{\pi }{W},k_y-n\frac{\pi }{W}\right)\cdot \sin {\mathrm{c}}_W(k_x,k_y)$$

$$Q_2=\sum _{m,n=-\infty }^{\infty }\left({(-1)}^m+{(-1)}^n\right)\delta \left(k_x-m\frac{\pi }{W},k_y-n\frac{\pi }{W}\right)\cdot \sin {\mathrm{c}}_W\left(k_x,k_y\right)$$

$$\sin {\mathrm{c}}_W(k_x,k_y)=\mathrm{sinc}(\frac{k_{x}W}{2\pi })\cdot \mathrm{sinc}(\frac{k_{y}W}{2\pi })$$

$$\mathrm{sinc}(k)=\frac{\sin (k\pi )}{k\pi }$$

$$\frac{{\tilde{E}}_k(u_x,u_y)\text{buildrel}.}{=\frac{A_{0}C}{\lambda f}(h_k{e}^{i\frac{2\pi }{\lambda f}{u}_y{y}_1}+\sum _{j=1}^{N_v}{W}_{jk}{v}_j/L\cdot e^{i\frac{2\pi }{\lambda f}{u}_y{y}_{j+1}}+b_k/L\cdot e^{i\frac{2\pi }{\lambda f}{u}_y{y}_{N_v+2}})\sin {\mathrm{c}}_W(0,0)\mathrm{sinc}\left(\frac{u_x{N}_{x}W}{\lambda f}\right)\mathrm{sinc}\left(\frac{u_y{N}_{y}W}{\lambda f}\right)}$$

$$I_k(0,0)=\frac{{(A_{0}C)}^2}{{(\lambda f)}^2}(h_k+b_k/L+\sum _{j=1}^{N_v}{W}_{jk}{v}_j/L{)}^2$$

$$I=\sum _k^{N_h}{I}_k(0,0)=\frac{{(A_{0}C)}^2}{{(\lambda f)}^2}\sum _{k=1}^{N_h}(h_k+b_k/L+\sum _{j=1}^{N_v}{W}_{jk}{v}_j/L{)}^2$$

$$\tilde{I}=N_h(N_v/L+1/L+1{)}^2\frac{I}{I_{norm}}=\sum _{k=1}^{N_h}(h_k+b_k/L+\sum _{j=1}^{N_v}{W}_{jk}{v}_j/L{)}^2$$

$$H_0=-L{\tilde{I}}_0=-2\sum _{i=1}^{N_h}\sum _{j=1}^{N_v}{v}_j{W}_{ji}-2\sum _{i=1}^{N_h}{b}_i-N_{h}L-\sum _{i=1}^{N_h}(\sum _{j=1}^{N_v}{v}_j{W}_{ji}{)}^2/L-\sum _{i=1}^{N_h}{(b_i)}^2/L-2\sum _{i=1}^{N_h}\sum _{j=1}^{N_v}{v}_j{W}_{ji}{b}_i/L$$

$$H_0=-2\sum _{i=1}^{N_h}\sum _{j=1}^{N_v}{v}_j{W}_{ji}-2\sum _{i=1}^{N_h}{b}_i-Q$$

$$H_k=-2\sum _{i\ne k}^{N_h}\sum _{j=1}^{N_v}{v}_j{W}_{ji}-2\sum _{i\ne k}^{N_h}{b}_i+2\sum _{j=1}^{N_v}{v}_j{W}_{jk}+2b_k-Q$$

$$p(\mathbf{v},\mathbf{h})=\frac{1}{Z}\exp [-H(\mathbf{v},\mathbf{h})/T],Z=\sum _{\mathbf{v},\mathbf{h}}\exp [-H(\mathbf{v},\mathbf{h})/T].$$

$$p(v_i|{\mathbf{v}}_{\backslash i},\mathbf{h})=p(v_i|\mathbf{h}),p(h_j|\mathbf{v},{\mathbf{h}}_{\backslash j})=p(h_j|\mathbf{v})$$

$$\begin{array}[]{ll}p({\mathbf{v}})&=\sum\limits_{\mathbf{h}}{p({\mathbf{v,h}})}=\frac{1}{Z}\sum\limits_{\mathbf{h}}{\exp[-H({\mathbf{v,h}})}]\\ &=\frac{1}{Z}\sum\limits_{\mathbf{h}}{\exp[\sum\limits_{ik}{{v_{i}}{W_{ik}}{h_{k}}+\sum\limits_{i}{{a_{i}}{v_{i}}}+}\sum\limits_{i}{{b_{i}}{h_{i}}}}]\\ &=\frac{{\exp({a^{T}})}}{Z}\sum\limits_{\mathbf{h}}{\exp[\sum\limits_{k}{{h_{k}}(\sum\limits_{i}{{v_{i}}{W_{ik}}}+}{b_{k}}})]\\ &=\frac{{\exp({a^{T}})}}{Z}\sum\limits_{\mathbf{h}}{\prod\limits_{k}{\exp[{h_{k}}(\sum\limits_{i}{{v_{i}}{W_{ik}}}+{b_{k}})]}}\\ &=\frac{{\exp({a^{T}})}}{Z}\sum\limits_{{h_{1}}}{\sum\limits_{{h_{2}}}{\cdots\sum\limits_{{h_{n}}}{\prod\limits_{k}{\exp[{h_{k}}(\sum\limits_{i}{{v_{i}}{W_{ik}}}+{b_{k}})]}}}}\\ &=\frac{{\exp({a^{T}})}}{Z}\prod\limits_{k}{\sum\limits_{{h_{k}}}{\exp[{h_{k}}(\sum\limits_{i}{{v_{i}}{W_{ik}}}+{b_{k}})]}}\\ \end{array}$$

$$p(h_j=1,\mathbf{v})=\frac{\exp (a^T)}{Z}\prod _{k\ne j}\{\sum _{h_k}\exp [h_k(\sum _i{v}_i{W}_{ik}+b_k)]\}\exp (\sum _i{v}_i{W}_{ij}+b_j)$$

$$p(h_j=1|\mathbf{v})=\frac{\exp (\sum _i{v}_i{W}_{ij}+b_j)}{\exp (\sum _i{v}_i{W}_{ij}+b_j)+\exp (-\sum _i{v}_i{W}_{ij}-b_j)}=\frac{1}{1+\exp (-2\sum _i{v}_i{W}_{ij}-2b_j)}$$

$$p(v_i=1|\mathbf{h})=\frac{1}{1+\exp (-2\sum _j{h}_j{W}_{ij}-2a_i)}$$

Full text

I Introduction

Restricted Boltzmann machine (RBM) is a stochastic generative neural network based on the Ising model, characterized by spin interactions solely between the visible and hidden layers [1, 2]. RBMs can be applied to a variety of tasks by learning the probability distribution of the training data. It is renowned for the development of deep learning [3, 4, 5, 6] and for demonstrating features in classification and recognition tasks, particularly generating new content [6, 7, 8, 9, 10]. Furthermore, by combining RBM with recurrent neural network, the time-varying RBM enables effective modeling of temporal data [11]. In RBMs, Gibbs sampling plays a crucial role by decomposing high-dimensional problems into multiple conditional distribution problems with Markov chains, and sampling sequentially from the conditional distribution [12, 13]. However, traditional electronic computing struggles with long chains of Gibbs sampling, especially when dealing with large-scale datasets, which demand enormous computational resources [6, 14, 15, 16].

Recently, there has been growing interest in developing unconventional computational architectures for simulating Ising Hamiltonians to accelerate specific computations. For example, to solve NP-hard problems, various Ising machines have been developed based on optical parametric oscillators [17, 18, 19, 20, 21], lasers [22, 23, 24, 25, 26], polariton [27, 28, 29], trapped ions [30], atomic and photonic condensates [31, 32], electronic memorisers [33], superconducting qubits [34, 35, 36], and nanophotonics circuits [37, 38, 39, 40, 41, 38, 42, 43, 44, 45]. In particular, the spatial photonic Ising machine (SPIM) has been demonstrated with reliable large-scale Ising spin systems, even up to ten-thousands of spins [46]. By leveraging multiple degrees of freedom of light, such as time [47], space [48, 49, 50], and wavelength [51], several division multiplexing schemes have been developed, harnessing the advantages of the ultrafast-speed and high-power-efficiency photonic computation [47, 51, 48, 49, 50]. Unlike traditional Von Neumann architectures, SPIMs encodes both the spin interaction matrix and spin configurations using spatial light modulators (SLMs) [46, 52, 53, 54, 47, 55, 56] or digital micromirror devices (DMDs) [57, 58], with measuring light intensity and instantaneously computing Ising Hamiltonian with Mattis-type interaction [46, 57, 52, 53], antiferromagnetic interaction [54], low-rank interaction [47, 55], four-body interaction [56]. However, all these Ising machines are designed to simulate Ising model wherein spins reside within a single layer and interact with each other, whereas RBMs partition spins into distinct visible and hidden layers, with spin interactions existing solely between two layers. Therefore, these computation schemes are unable to accelerate the computation for RBMs, due to the discrepancy of the spin distribution and interaction.

In this article, we propose a photonic restricted Boltzmann machine (PRBM) that leverages photonic computing to accelerate Gibbs sampling and demonstrate content generation tasks. The PRBM is designed based on wavelength-division multiplexing spatial Ising machines, inheriting the scalability of SPIM architecture. By developing an encoding method, we eliminate the need for the decomposition computation of interaction matrix and reduce the computational complexity of Gibbs sampling from $O(N)$ to $O(1)$, thereby significantly reducing computational costs. In addition, with the advantage of photonic computation of non-Von Neumann computing architecture, we circumvent the memory storage of interaction matrix, substantially enhancing computation efficiency for large-scale RBM. We further validate Gibbs sampling process by simulating analytically solvable two-dimensional Ising model [59] and experimentally observe the phase transition at the critical temperature agreeing well with exact theoretical value. Moreover, we experimentally demonstrate that the proposed PRBM effectively acquires the capability of generating new content, with successful implementations in image generation and temporal data processing, as exemplified by music generation. Finally, we prospect the scalability and the training advances of the PRBM framework, suggesting its potential to offer a promising pathway for applications of generative artificial intelligence (AI).

II Photonic computing scheme of restricted Boltzmann machine

We consider RBMs with a Hamiltonian $H({\bf{v,h}})=-\sum\limits_{ij}{{v_{i}}{W_{ij}}{h_{j}}-\sum\limits_{i}{{a_{i}}{v_{i}}}-}\sum\limits_{i}{{b_{i}}{h_{i}}}$, where ${v_{i}}$ and ${h_{i}}$ taking values of $\pm 1$ represent the spin configurations in the visible and the hidden layers. The spins interact between the two layers solely as $W_{ij}$, where ${a_{i}}$ and ${b_{i}}$ are the magnetic field terms for visible layer spins and hidden layer spins, respectively (see Figure 1(a) and the Methods). We note that SPIMs simulating such spin interactions in RBMs has to introduce substantial redundant zero-value elements [60] and further requires digital computation of eigenvalue decomposition [47, 48] or Cholesky-like decomposition [51] with the cost of computation complexity $O({N^{3}})$. Such computation cost becomes significant when the total spin number in the visible and hidden layers is large, and is further intensified when accounting for external magnetic fields on each spin. Additionally, SPIMs utilize Metropolis-Hastings (MH) sampling evaluating the impact of a single spin flip on the global system, whereas Gibbs sampling evaluates the impact of the locally correlated system on a single spin, which is crucial for accelerating computation in complex high-dimensional problems. Therefore, how to effectively encode the visible and hidden layers and how to implement Gibbs sampling constitute two key challenges in demonstrating the photonic computing acceleration for RBMs.

Figure 1(b) shows the proposed PRBM by the photonic computing of Gibbs sampling, following the measurement and feedback scheme. Here a collimated supercontinuum laser is diffracted by a grating and different wavelength light focuses on a SLM along the $x$ axis with a cylindrical lens. The modulated light beam reflected from SLM is transformed by a lens, resulting in an incoherent intensity summation of different wavelengths and coherent interference for each wavelength at the back focal plane (detailed information is provided in Supplementary Materials Section $\rm I$ ).

To enable the photonic computing of Gibbs sampling, we divide the SLM into three distinct regions along the $y$ axis. We take the example of the generation process as illustrated in Figure 1(c). When sampling from the visible layer ${\mathbf{v}}$ to the hidden layer ${\mathbf{h}}$, for each hidden-layer spin ${h_{k}}$, the visible-layer spin ${{v_{i}}}$ with the spin interaction ${W_{ik}}$ and under the magnetic field ${b_{k}}$ are encoded at the $k$-th wavelength, and the phase modulation along the $y$ axis of the SLM are given as

[TABLE]

Here an auxiliary spin $s=1$ is used for magnetic field ${b_{k}}$ in Region ${\rm III}$. We encode each spin by a macro-pixel which consists of ${N_{x}}\times{N_{y}}$ pixels, where $m$ and $n$ are the pixel indices within each spin, $1\leqslant m\leqslant{N_{x}}$ and $1\leqslant n\leqslant{N_{y}}$.

In order to implement the spin interaction and the external magnetic field, we propose the Gauge transform [52, 54]: In Region ${\rm II}$ and ${\rm III}$, the phase modulations have additional contributions from a checkerboard modulation, corresponding to rotating each spin by the angle ${\alpha_{ik}}=\arccos({W_{ik}}/L)$ and ${\beta_{k}}=\arccos({b_{k}}/L)$, which are equivalent spins $v_{ik}^{{}^{\prime}}={v_{i}}{W_{ik}}$ and $s_{k}^{{}^{\prime}}=s{b_{k}}$, respectively. The normalization parameter $L$ is the maximum absolute value of the interactions and external magnetic fields.

With this encoding method, we perform Gibbs sampling in two steps. In the first step, we set all hidden layer spins as ${h_{i}}=1$ and encode them with visible layer spins onto the SLM following Equation 1. Thus, the normalized intensity ${\tilde{I}_{0}}$ at the center position of the back focal plane is obtained, and the Hamiltonian is computed as ${H_{0}}=-L{\tilde{I}_{0}}=-2\sum\limits_{i=1}^{{N_{h}}}{\sum\limits_{j=1}^{{N_{v}}}{{v_{j}}{W_{ji}}-}}2\sum\limits_{i=1}^{{N_{h}}}{{b_{i}}-}Q$, where $Q$ is a constant, ${N_{v}}$ and ${N_{h}}$ are the number of the spins in the visible and hidden layers, respectively (Supplementary Materials Section $\rm{II}$). In the second step, every spin in the hidden layer is flipped independently, as indicated by the gray shadowed box in Figure 1(c). For the $k$-th spin, we set ${h_{k}}=-1$ and measure the intensity corresponding to the Hamiltonian ${H_{k}}=-2\sum\limits_{i\neq k}^{{N_{h}}}{\sum\limits_{j=1}^{{N_{v}}}{{v_{j}}{W_{ji}}-}}2\sum\limits_{i\neq k}^{{N_{h}}}{{b_{i}}+2\sum\limits_{j=1}^{{N_{v}}}{{v_{j}}{W_{jk}}}+2{b_{k}}-Q}$. Therefore, by measuring the difference of intensity in two steps, we can compute $\Delta H_{k}$ as

[TABLE]

Based on the feedback scheme, we update the spin configuration with Gibbs sampling at temperature $T$, by calculating the probability of the $k$-th spin being $+1$ as $P({h_{k}}=1|{\bf{v}})=1/[1+\exp(\Delta H_{k}/T)]$.

We note that the proposed encoding method for Gibbs sampling does not require the decomposition of the spin interaction matrix [47, 48, 51] and can be generalized for sparse spin interactions (see more details in the Methods). When sampling from the hidden layer to the visible layer, we encode the visible layer spin ${v_{k}}$ and the external magnetic field ${a_{k}}$ in region ${\rm I}$ and ${\rm III}$, respectively, and the equivalent spin $h_{ik}^{{}^{\prime}}={h_{i}}{W_{ki}}$ in Region ${\rm II}$. In the process of generating an image, we begin with a randomized visible layer spin configuration encoded onto Region II. By performing iterative photonic computing of Gibbs sampling between the visible and hidden layers, a generated content is decoded from the final visible layer as shown in Figure 1(c).

III Validation of PRBM by observing phase transition

We first validate the proposed PRBM with Gibbs sampling by observing the phase transition in the two-dimensional Ising model $H=-J\sum\limits_{<ij>}\sigma_{i}\sigma_{j}$, where each spin interacts only with its four nearest neighbors of the strength $J$. Analogous to RBM, we transform the spins at the four corner positions of the square into visible layer (red circular dots in Fig. 2(a)) and hidden layer (blue triangles in Fig. 2(a)), respectively. Through this transformation, we study the phase transition and compare the experimental results with the analytical critical temperature [59]. We experimentally conduct photonic computing of the Gibbs sampling on a $10\times 10$ lattice at fourteen different temperatures. In order to calculate the phase diagram, we generate 20 samples at each temperature, with each sample obtained after 10 iterations from a random initial spin configuration, with each complete iteration consisting of Gibbs sampling from the visible layer to obtain the hidden layer, and vice versa. Figure 2(b) displays the obtained spin configurations after 10 iterations at four different temperatures $T=6.1J,3.1J,1.9J$, and $1.1J$. At high temperature, such as $T=6.1J$, the spin configuration exhibits disordered, while as the temperature decreases, cluster of aligned spins begin to emerge and the spin configuration becomes increasingly ordered.

To precisely evaluate the phase transition temperature, we calculate the heat capacity $C=\frac{{\left\langle{{H^{2}}}\right\rangle-{{\left\langle H\right\rangle}^{2}}}}{{N{T^{2}}}}$ and susceptibility $\chi=\sum\limits_{i=1}^{N}{\frac{{1-m_{i}^{2}}}{{NT}}}$, where ${m_{i}}$ represents the statistical average of a spin and $N$ is the number of spins. As shown in Figure 2(c), the experimental results show that the phase transition temperature is ${T_{c}}=2.3J$, which agrees well with the theoretical predictions ${T_{c}}=\frac{{2J}}{{\ln(1+\sqrt{2})}}\approx 2.27J$ in the infinite 2D model [59]. We also employ the annealing method [61] with Gibbs sampling to search the ground state (Supplementary Materials Section $\rm III$). All these results demonstrate that, despite limited noise and aberrations, the proposed photonic computing of Gibbs sampling effectively simulates the spin interaction between the visible and hidden layers for RBMs.

IV Image content generation and restoration by PRBM

Next we demonstrate the capability of the photonic restricted Boltzmann machine for content generation. In image generation, we utilize the Fashion MNISIT dataset [62] and MNIST digit images [63] as the visible layer spins to train the RBM, where each pixel in an image corresponds to a spin. Before training, we compress the original images with $28\times 28$ spins to a resolution of $14\times 14$ (yielding 196 spins) and binarize the images. During both the training and generation processes, Gibbs sampling is conducted at a fixed temperature $T=1$. Furthermore, we train the model in the absence of magnetic field and separately train it on the fashion items ‘Boot’, ‘Pants’, and the digit ‘0’. Following the Hinton recipe [64], we initialize the interaction weights with small random values drawn from a zero-mean Gaussian distribution (Supplementary Materials Section $\rm{IV}$).

To generate an image by PRBM, we first randomly initialize the state vector ${{\mathbf{v}}^{0}}$ and subsequently obtain ${{\mathbf{v}}^{K}}$ as the terminal generated image after $K=15$ iterations by performing photonic computing of Gibbs sampling. With the stabilized spin configuration in the visible layer, Figure 3 illustrates three distinct image generation processes for each experimental item. Additional experimental results on image generation are provided in Supplementary Materials Section $\rm{VII}$. The diversity observed across these generated images demonstrates the PRBM’s capacity to produce new content with significant variability.

We further validate photonic computing of Gibbs sampling by evaluating whether the trained model exhibits overfitting. Here we use PRBM to perform restoring-image test, where the initial images are excluded in the training dataset and distorted. We compress images of the fashion item ‘Boot’ and the digit ‘0’ to $20\times 20$ spins and train the model. For the image restoration experiments, we mask three different ‘Boot’ images randomly selected from the test dataset, as indicated by the red areas in Figure 4(a). Figure 4(b) represents the initial input images ${{\mathbf{v}}^{0}}$ (the first column), where spins in the masked regions are set to a value of $-1$. After $K=15$ iterations, the three restored images ${{\mathbf{v}}^{K}}$ are obtained, as shown in the last column of Figure 4(b).

As another example, we also restore three different ‘0’ images from the test dataset by adding random noise, as indicated by the red marks in Figure 4(c). Before performing image restoration, we set the spins in the red-marked regions to the negative values of their respective original values, as is shown in the first column ${{\mathbf{v}}^{0}}$ of Figure 4(d). After 15 iterations, the final restored images ${{\mathbf{v}}^{K}}$ are presented in the last column of Figure 4(d). More other experimental results about image restoration are shown in Supplementary Materials Section $\rm{VII}$). All these results demonstrate that the image generation capabilities of the PRBM are not attributable to overfitting and also validate the proposed photonic computing of Gibbs sampling.

V Temporal content generation by PRBM

To demonstrate time-varying RBM systems, we generate pieces of piano music by using Recurrent Neural Network-Restricted Boltzmann Machine (RNN-RBM) [11], which forms a sequential step through each RBM cell (the green box in Figure 5(a)) along the temporal dimension. Different to the ordinary RBM, the magnetic field coefficients in RNN-RBM are dynamically updated at each time step, whereas the interactions parameters remain invariant. In the experimental implementation of piano music generation with the PRBM, as the time step $t$ increases, we only need to modify the rotation ${\beta_{k}}$ in Region ${\rm III}$ of the SLM by utilizing the Gauge transform. We use the piano dataset Nottingham [65] as the training dataset. Within the PRBM framework, a total of 88 piano keys ranging from A0 to C8 are mapped onto 88 spins $v_{i}$ in the visible layer of RBM, where $v_{i}$ assumes the value of +1 to represent the pressing of the $i$-th piano key and $-1$ otherwise.

Figure 5(a) demonstrates the music generation process. The gray arrows indicate the connections between the elements. The red arrows signify that the spins ${\bf{v}}_{t}^{K}$ in each RBM cell correspond to the specific piano keys requiring activation at the current time step. In the experiments, the proposed PRBM executes the music generation composed of 150 time steps, where each time step stands for a quaver. The architecture of each RBM cell within this temporal sequence incorporates 88 spins in the visible layer, corresponding to the full range of piano keys, and 96 spins in the hidden layer to capture the underlying musical patterns and dependencies.

Figure 5(b) presents one of the pieces in training dataset. Following the training process, the rhythmic patterns and stylistic characteristics of the music are encoded within the coefficient parameters. Subsequently, we generate new melodies exhibiting similar rhythmic and stylistic properties. For each time step $t$, ${\bf{v}}_{t}^{0}$ is randomly initialized as the spins in visible layer of RBM. The generated music ${\bf{v}}_{t}^{K}$ is produced after 20 iterations, as shown in the $t$-th column in the Fig. 5(c). The ${\bf{v}}_{t}^{K}$ is then utilized for updating the RNN elements ${{\mathbf{r}}^{t}}$ in current time step. The generation processes are conducted repeatedly, and after 150 time steps, a piano music piece is generated as shown in Fig.5(c). The generated musical composition demonstrates remarkable similarity to the training examples in terms of rhythmic structure and stylistic attributes, thereby validating the efficacy of the proposed PRBM for temporal content generation.

VI Discussion and Outlook

In conclusion, we propose a photonic restricted Boltzmann machine with introducing an encoding method to accelerate Gibbs sampling processes. In contrast to SPIM, the proposed encoding method eliminates the decomposition computation of interaction matrix and magnetic field, thereby significantly improving computing efficiency. We verify the PRBM and the photonic computing of Gibbs sampling by observing the phase transition in the 2D lattice lsing model. Furthermore, we experimentally demonstrate the capability of PRBM to generate new content by generating and restoring images. To further illustrate the advantages of the encoding method, we demonstrate that PRBM can implement time-varying spin interactions for temporal content generation, such as music, where non-decomposition computation is particularly crucial for efficient processing.

We note that the proposed system effectively accelerates computation through parallel processing and reduces the computational complexity of Gibbs sampling from $O(N)$ in the traditional electronic computing to $O(1)$ in the photonic computing. Moreover, with the non-Von Neumann photonic computing architecture, the interaction matrix and the magnetic fields are encoded on the SLM, significantly saving the memory cost and breaking the data transportation bottleneck between CPU/GPU and memory in the traditional electronic computation, especially substantially enhancing computation efficiency for large-scale RBM.

We also evaluate the computing capacity and scalability of the proposed system. For each photonic computing of Gibbs sampling, the computational load is represented by Equation (2), which involves $N$ real-number multiplication and $N$ summation, where $N$ is the number of spins. The time per step (TPS) in the present experiment is primarily limited by the modulation speed of the liquid-crystal SLM. With the development of various cutting-edge SLMs capable of gigahertz modulation rate and low power consumption[66, 67, 68], the minimum TPS could scale down to ${10^{-9}}$s, as PRBM requires only a single pixel photodetector. Moreover, the PRBM can be scaled to larger spin size by extending the wavelength range and increasing the number of SLM pixels. As the pixel scale of the modulator reaches the wavelength scale [68], it is possible to achieve $10^{10}$ pixels in a $100\rm{mm}\times 100\rm{mm}$ area, enabling the realization of models with 10 billion parameters. In this case, the number of spins is $N=10^{5}$, the proposed experimental framework can achieve a maximum of 200-tera-floating-point operations per second (TFLOPS).

It is worth noting that using the proposed photonic Gibbs sampling method can potentially improve computing efficiency and reduce resource consumption for training processes. Current generative AI requires tremendous computation resources for training. For example, Transformer-based language models using digital computers require computational operations on the order of $6MBS$, where $M$, $B$ and $S$ are the number of model parameters, the training batch size and the iteration number, respectively [69]. In comparison, the total computation operations using digital computers for RBM are the same order of magnitude as that of large language models, approximately $9MBS$, where $M\approx{N^{2}}/4$, $N$ is the total spin number (Supplementary Materials Section $\rm{V}$). Owing to the accelerated Gibbs sampling enabled by the proposed PRBM, the complexity can be reduced to $3MBS+3NBS$. Therefore, photonic computing acceleration enables the training efficiency of RBM outperforms the Transformer-based models, which is a significant cue to develop unconventional computational architectures for generative AI.

With the spin size $N$ increasing, the operations completed per second grow, which indicates that the computational capacity and training efficiency of the proposed PRBM could outperform electronic computing. Supplementary Materials Table S1 shows that PRBM reduces the training time by about two orders of magnitude compared to the NVIDIA H100 when training the GPT-3 model with the same number of parameters and training data. Given that the proposed photonic computing scheme demonstrates strong compatibility with diverse data types, including both sequential and non-sequential data, we believe that PRBM could achieve high practical efficiency for learning complex probability distributions and generating sophisticated contexts, including potential applications in language models.

Methods

Restricted Boltzmann machine At temperature $T$, the probability of a spin configuration follows the Boltzmann-Gibbs distribution $P({\bf{v,h}})\propto\exp[-H({\bf{v,h}})/T]$, where $H$ is the energy of that state (here we take the convention of the Boltzmann constant ${k_{B}}=1$). As shown in Figure 1(a), the training set is binarized as the spins in the visible layer, and the training objective of RBMs is to maximize the log-likelihood function $\ell({\bm{\theta}}){\rm{=}}{\left\langle{\ln P({\bf{v,h}})}\right\rangle_{data}}$ with the parameters ${\bm{\theta}}=\left\{{{W_{ij}},{a_{i}},{b_{i}}}\right\}$, at temperature $T=1$. Gibbs sampling is utilized to compute the partial derivatives of $\ell({\bm{\theta}})$ with respect to each parameter, which enables the update of parameters. The probability of the $k$-th spin being $+1$ in the hidden layer is $P({h_{k}}=1|{\bf{v}})=P({h_{k}}=1,{\bf{v}})/P({\bf{v}})=1/[1+\exp(-2({b_{k}}+\sum\nolimits_{i}{{v_{i}}{W_{ik}}}))]$ (Supplementary Materials Section $\rm{VI}$), while the probability of the $k$-th spin in the visible layer being $+1$ is calculated as $P({v_{k}}=1|{\bf{h}})=1/[1+\exp(-2({a_{k}}+\sum\nolimits_{i}{{W_{ki}}{h_{i}}}))]$. In the generation process, Gibbs sampling is employed to generate new content of the same type based on the trained parameters.

General sparse interaction matrix The proposed encoding method is also applicable to Ising models with complex sparse interactions. During Gibbs sampling, the spin ${\sigma_{i}}$ is encoded in Region ${\rm I}$ of SLM, while the spins connected to ${\sigma_{i}}$ are encoded in Region ${\rm II}$. Therefore, this photonic computing system can be extended to simulate general Ising models with sparse interactions. Especially for combinatorial optimization problems, such as Graph Partitioning Problem with arbitrary sparsity, the system is expected to exhibit unique performance without requiring matrix decomposition.

Acknowledgement

The authors acknowledge funding through the National Natural Science Foundation of China (124B2085, 12474393 and 12174340), and the National Key Research and Development Program of China (2022YFA1405200).

I Experimental set-up for photonic restricted Boltzmann machine

In Figure S1, a supercontinuum laser source (Anyang SC-5) produces a collimated Gaussian beam. The beam waist radius undergoes tenfold expansion through a lens system comprising L1 (50 mm in focal length) and L2 (500 mm in focal length). The expanded beam subsequently illuminates a reflective diffraction grating (inscribed line density of 600/mm), and the cylindrical lens (150 mm in focal length) focus the light with different wavelengths onto an SLM along the $x$-axis, while the $y$-axis pixels are illuminated coherently by the same wavelength. Polarizer P ensures linear polarization alignment with the SLM’s long display axis. L3 (100 mm in focal length) and L4 (200 mm in focal length) magnify twofold at the SLM plane (Holoeye PLUTONIR-011). Lens FL (200 mm in focal length) performs a Fourier transformation of the optical field, resulting in an incoherent intensity summation of different wavelengths and coherent interference for each wavelength at the back focus plane. An sCMOS camera (Andor Zyla 5.5) positioned at this focal plane records the intensity distribution, following the measurement and feedback scheme. To account for finite SLM pixel dimensions and sensor resolution limitations, intensity measurements are obtained through spatial integration around central regions rather than single-point detection. In the experiment, the operational wavelengths are selected from 588 $\rm nm$ to 611 $\rm nm$. Each individual wavelength projects a beam covering approximately two pixels along the $x$-axis of the SLM, approximately 13.5 $\rm\mu m$. Each pixel of SLM has a physical width of 8 $\rm\mu m$. Based on the grating equation, the wavelength interval for two neighboring wavelengths is estimated about $\Delta\lambda=0.053\rm{nm}$.

II Photonic computing of Gibbs sampling

We take the example of the generating process as illustrated in Figure 1(c) in the main text. We first assume that the light intensity for each wavelength on the SLM is uniform and equal to $A_{0}$. When sampling from the visible layer to the hidden layer, for each wavelength, the phase modulation $\phi_{j,k}^{m,n}$ on the SLM is defined as Eq. (1) in the main text. With proposed Gauge transform, each spin encoded by ${N_{x}}\times{N_{y}}$ pixels. Therefore, the electric field modulated after SLM is expressed as

[TABLE]

where ${G_{+}}$ and ${G_{-}}$ are written as

[TABLE]

In this context, $(x,y)$ represent the spatial coordinate on the SLM, and $y_{j}=j{{N}_{y}}W$ represents the position of the $j$-th spin in $y$ direction on SLM. The rectangular function is defined as

[TABLE]

where $W$ represents the pixel width of the SLM. $H$ is the checkerboard function, so $1-H$ is the complementary function. The checkerboard function is defined as

[TABLE]

where $m$ and $n$ are two integers, and the mathematical symbol $\otimes$ denotes the convolution operation. The optical field at the focal plane after passing through the lens is expressed as the result of a two-dimensional Fourier transform, written as follows

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Due to the checkerboard modulation, the optical beam is diffracted into multiple orders at $(m\frac{\pi}{W},n\frac{\pi}{W})$ in the angular spectral space. Furthermore, we convert the angular spectrum coordinates $({u_{x}},{u_{y}})$ of the detection area with ${u_{x}}={k_{x}}\frac{{f\lambda}}{{2\pi}}$ and ${u_{y}}={k_{y}}\frac{{f\lambda}}{{2\pi}}$. Due to the minimal overlap of the light fields between different orders, it can be neglected. The zeroth-order diffraction field, where both $m$ and $n$ are equal to zero, can be expressed as

[TABLE]

where $C={N}_{x}W\cdot{{N}_{y}}W$. Therefore, the detected intensity at ${u_{x}}=0$ and ${u_{y}}=0$ is

[TABLE]

Owing to the incoherence among different wavelengths, the overall detection intensity at the zeroth order is given by

[TABLE]

where $N_{h}$ is the number of the spins in the hidden layer. We define all coefficients and spins to be equal to 1, with detected intensity ${I_{norm}}=\frac{{{{({A_{0}}C)}^{2}}}}{{{{(\lambda f)}^{2}}}}{N_{h}}{({N_{v}}/L+1/L+1)^{2}}$. By normalizing the detected intensity to $I_{norm}$, the normalized intensity $\tilde{I}$ is written as

[TABLE]

Then we take the example of the generating process, we perform photonic computing of Gibbs sampling in two steps. For example, when sampling from the visible layer to obtain the hidden layer, in the first step, we set all hidden layer spins ${h_{i}}=1$ and encode them with a random initial state onto the SLM following Equation (1) in the main text. According to Eq.(S8) and Eq.(S9), the normalized intensity ${\tilde{I}_{0}}$ at the center position of the back focus plane is obtained as ${\tilde{I}_{0}}=\sum\limits_{k=1}^{{N_{h}}}{{\rm{(1+}}{b_{k}}{\rm{/}}L{\rm{+}}\sum\limits_{j=1}^{{N_{v}}}{{W_{jk}}{v_{j}}{\rm{/}}L{{\rm{)}}^{2}}}}$. Thus, the Hamiltonian is computed as

[TABLE]

When calculating the Gibbs sampling probability for the hidden layer spins, the visible layer spins and the magnetic field remain unchanged. Thus, the above expression simplifies to

[TABLE]

where $Q$ is a constant. In the second step, each spin in the hidden layer is flipped independently. For the $k$-th spin, we set ${h_{k}}=-1$ and compute the corresponding Hamiltonian as

[TABLE]

Thus, the Hamiltonian difference is $\Delta H_{k}=({H_{0}}-{H_{k}})/2=-2\sum\limits_{j=1}^{{N_{v}}}{{v_{j}}{W_{jk}}}-2{b_{k}}$, and the Gibbs sampling probability is $P({h_{k}}=1|{\bf{v}})=1/(1+{{\mathop{\rm e}\nolimits}^{\Delta H_{k}}})$. The photonic computing system we propose reduces the computational complexity of Gibbs sampling from $O(N)$ to $O(1)$, substantially enhancing computation efficiency for large-scale RBM.

III Ground-state search for two-dimensional Ising model

Figure S2 illustrates the variation of magnetization $M$ (orange line) and Hamiltonian $H$ (blue line) during the cooling process from high to low temperature. The black line denotes the temperature variation during the iterative process, which undergoes a decay from an initial temperature $T=6.1J$ to a final temperature $T=0.3J$. The Hamiltonian gradually decreases with decreasing temperature, eventually reaching its minimum value. During the cooling process, the magnetization shows that the spin configuration transitions from a paramagnetic phase to a ferromagnetic phase.

IV Iterative generating with photonic computing of Gibbs sampling

In the image generation experiments, we shrink the original dataset to $14\times 14$ pixels and binarize them as the visible layer spins. We set the number of hidden layer spins as 196. To train the model, we use the contrastive divergence (CD-1) algorithm and train the model in the absence of magnetic field. We set the learning rate for the interaction parameters to 0.003, and to avoid overfitting, the weight decay parameter is set to 0.0002. The initial values of the interaction parameters ${W_{ij}}$ follow a Gaussian distribution with a mean of 0 and a variance of 0.01. We train on 1000 images, dividing them into 50 batches for training. The total number of training epochs is 2000. During the training process, we introduce the Momentum method to address issues such as slow convergence, oscillations, and getting trapped in local optima in the gradient descent algorithm. Figure S3(a) shows the evolution of the log-likelihood function in the training process for item ‘Boot’. The log-likelihood function gradually increases and approaches 0 as the number of training iterations increases, demonstrating that our training is effective. In the generation of the image, an initial state is randomly generated, and Gibbs sampling is performed for each spin of the hidden and visible layers in turn, and so on for 15 iterations, with the final visible layer spins representing the generated image.

In the experiments of image restoration, we shrink the original dataset to $20\times 20$ pixels and binarize them as the visible layer spins. The number of hidden layer spins is same to that of visible layer spins. In the process of learning the parameters, we employ the same learning rate and initial values for the parameters as used in the generation experiments. We train on 1500 images, dividing them into 50 batches for training. The total number of training epochs is 2000. Figure S3(b) shows the evolution of the log-likelihood function in the training process for item ‘Boot’. To test over-fitting, we use different damaged validation images as the initial states and as the spins of visible layer. The Gibbs sampling is performed for each spin of the hidden and visible layers in turn, and so on for 15 iterations, with the final visible layer spins representing the restored images.

For temporal content generation in the PRBM, we use the piano dataset (Nottingham) as the training dataset. The RNN-RBM model is trained with CD-15 algorithm to maximize the log-likelihood function, and we introduce the Adam optimization algorithm during the training process. All parameters are initialized according to a Gaussian distribution with a mean of 0 and a variance of 0.01, and the learning rate is set to 0.001. We train on 1000 piano songs, dividing them into 10 batches for training. The total number of training epochs is 200, so the training iteration is 2000. Figure S3(c) illustrates the evolution of the log-likelihood function.

V Derivation of the digital computational complexity for training RBM

Training an RBM typically uses the CD-$k$ (Contrastive Divergence) algorithm, where in most cases $k=1$. Here, we also use $k=1$ to calculate the computational complexity. Assuming both the hidden layer and visible layer have $N/2$ spins, and the number of training parameters including interactions and magnetic field is $M$, where $M={N^{2}}/4+N$. Since the number of magnetic field is negligible compared to the number of interactions, we define $M\approx{N^{2}}/4$. For given visible layer, calculating the probability of each spin in the hidden layer using Gibbs sampling requires $N/2$ real-number multiplication, $N/2$ summation and one computation of the sigmoid function. Since the computational cost of the sigmoid function is very small compared to the multiplication and summation, it can be neglected. Therefore, the computational complexity for calculating all the spins in the hidden layer is $N\times N/2$. Similarly, the computational complexity for reconstructing the visible layer from the hidden layer is $N\times N/2$. Based on CD-1 algorithm, an additional computation is required to sample the hidden layer from the reconstructed visible layer, which also has a computational complexity of $N\times N/2$. Assuming that the batch size for each parameter update is $B$, the total digital computational complexity required for Gibbs sampling is the summation of the aforementioned calculations, which yields a final complexity of $3{N^{2}}BS/2$ with total training iteration $S$.

When calculating the gradient of the interactions for batch size $B$, $2MB$ multiplication, and $MB$ summation operations are required in each iteration. The computation of the magnetic field gradient requires only $2NB$ summation, which is negligible compared to the computational cost of the interactions. Therefore, the total computation complexity required to train the RBM is approximately $C\approx 9MBS$, and Gibbs sampling constitutes a significant portion of the total computational complexity. With the proposed PRBM, the computational complexity of Gibbs sampling is reduced to $3NBS$. Therefore, the total computational complexity for training the RBM is $3MBS+3NBS$.

VI Derivation of Gibbs sampling

In general, Gibbs sampling involves the visible and hidden layer spins taking values of 0 or 1. However, given that the spin in Ising model is -1 or 1, the original formula for calculating the Gibbs sampling probability cannot be applied to the system proposed in this article. Therefore, we present the derivations of Gibbs sampling in the proposed system.

Given the Hamiltonian $H$, the joint probability distribution of the model is represented as

[TABLE]

Furthermore, we set the temperature $T$ in the experiment at 1. In the process of Gibbs sampling, the conditional probability of each spin is given by

[TABLE]

It is known from the conditional probability formula that $p({h_{j}}=1|{\mathbf{v}})=\frac{{p({h_{j}}=1,{\mathbf{v}})}}{{p({\mathbf{v}})}}$. The marginal probability can be expressed as

[TABLE]

So

[TABLE]

Therefore

[TABLE]

Similarly

[TABLE]

VII Experimental results on image generation and restoration

We have trained the RBM on all digit classes from the MNIST dataset, as well as the ‘Boot’, ‘Pants’, and ‘T-shirt’ categories from the Fashion-MNIST dataset. The experimental generation results for digit ‘0’, ‘Boot’, and ‘Pants’, obtained through our proposed PRBM and encoding method, are presented in the main text. The following figures demonstrate the experimental generation processes for the remaining image categories. In Figures S4-S13, the first column displays randomly initialized spin configuration, the subsequent columns (from left to right) illustrate the iterative stages of the image generation process, and the final column presents the generated image. Figure S14 demonstrates the image restoration process with three different images, which are sourced from the test set. The first column displays the corrupted image with a masked region, which serves as input for the restoration experiments. Similarly, Figure S15 demonstrates the restoration process for three images corrupted by random noise.