1. Introduction

The field of optical computing is witnessing a rapid development thanks to its unique properties such as high parallelization, low latency and essentially low energy cost [1]. Within various optical computing configurations [2–4], the optical linear transform (OLT) technique plays a crucial and fundamental role. The evolution of optical complex amplitudes in these systems is linear in nature and can be formulated using a complex-valued transformation matrix (TM) [5]. The OLT has been realized both in free space and on waveguide-based platforms. In free space, implementations based on meshes of beam splitters and phase shifters [6,7], cascaded phase masks [8] or diffractive layers [9] have been pushed forward. For implementing the OLT in waveguides, programmable Mach–Zehnder interferometer (MZI) [10], multimode interference (MMI) coupler [11] or multiport directional coupler [12] meshes have been proposed. In addition to these carefully fabricated and/or aligned optical architectures, the OLT has also been implemented using complex media, such as chaotic cavity [13–15], diffuser in free space [16] and waveguides [17,18]. Basically, the design of a user defined OLT system adopts two methodologies. One is the deep learning-based method which requires a pre-known database for training, validations and tests [9]. This method spends a substantial amount of computing resources preparing the learning data, as well as updating a large number of weight values in each iteration. The other methodology employs the general-purpose optimization algorithms to iteratively search for solutions [10–12,16]. Although these methods do not require learning from a database, the efficiency still leaves much room for improvement as the program does not necessarily follows the steepest converging path. In general, the current schemes for designing the required OLT need a highly efficient algorithm specially designed for obtaining the configurations of optical components from a target TM.

On the other hand, one basic and important issue in optical computing is to create massive optical complex data as the very input to the large-scale, parallel optical processing networks [19]. Previously, various methods have been proposed for encoding the complex optical field onto phase-only spatial light modulator (SLM) [20–23]. However, these techniques have two major drawbacks. Firstly, when generating multiple complex beams using phase gratings [24], controlling the amount of light diffracted into one direction will inevitably influence the beams in other positions. Hence, the beams in the optical array are not completely independent. Secondly, the advanced phase profile generation methods often require specific optical setups including Fourier lenses, various kinds of spatial filters [25,26], or multiple phase masks [27]. Therefore, a straightforward way based on the OLT design to generate optical complex data/beams with minimal number of physical components is highly desired.

In this work, we implement the OLT system using a single or multiple cascaded phase mask(s). A learning-free iterative algorithm is proposed to obtain the phase masks from a given complex-valued linear matrix. By comparing with the widely adopted deep learning-based methods and general-purpose algorithms, the proposed algorithm shows its superiority by eliminating redundant mathematical operations. With this new algorithm, we apply the OLT construction technique to generate arrays of optical complex beams. The experiments verify that both the amplitude and phase of each beam in the array can be separately controlled without affecting the others. It is also demonstrated that the optical system can be as simple as a light source and a single SLM without multiple reflections. Neither physical spatial filter nor Fourier lens is required. Both numerical calculations and experiments confirm that the number of required phase masks can be reduced to one, so long as it provides sufficient total pixels. Our method ensures the parallel generation of massive optical complex data using a simple optical configuration and a highly efficient algorithm. Both merits are attractive for optical information processing applications, such as input gate for optical computers [28,29] and generator for high-dimensional quantum optical states [30,31].

2. Design of the optical linear transform

To begin with, we assume that the two-dimensional (2D) input optical field I and output optical field O have the dimensions of n × n and m × m, respectively. Both I and O are composed of the point sampled optical complex amplitudes. In order to link the complex-valued data of these two optical fields with a 2D matrix, they are first reshaped into two vectors $\vec{X}$ and $\vec{Y}$ with the dimensions of n ²  × 1 and m ²  × 1, respectively. Here, the element in the vectors, $\vec{X}(g )$ or $\vec{Y}(g )$ is associated with the pixel of the optical fields I(a, b) or O(a, b) through g = (a − 1)n + b or g = (a − 1)m + b. Afterwards, the relation between the input and output fields can be represented as:

In Eq. (2), T ₀ and T _N are the TM representing the free-space propagation from the input beam to the first SLM and that from the last SLM to the output plane, respectively. T _i (i = 1 … N−1) is the complex matrix describing the free-space propagation from the i-th plane to the (i+1)-th one. These TMs can be derived from the Rayleigh-Sommerfeld diffraction formula [33,34]. In the matrices presenting the SLMs, φ_{i ,j} is the phase delay imposed by the j-th element on the i-th phase mask, and k is the total number of phase elements on each phase mask.

 

Fig. 1. The matrix representation of a general OLT system based on cascaded phase masks. The input and output optical complex fields of the OLT system I and O are complex-valued matrices representing the arrays of sampling points of the optical complex amplitudes. x_ij and y_ij are the complex amplitudes of the optical field at each sampling point in the input and output planes, respectively. The matrices I and O are reshaped into two complex-valued vectors $\vec{X}$ and $\vec{Y}$, respectively. The desired OLT process can be realized by designing the phase profiles on the SLMs.

Download Full Size | PDF

Previously, the optical configuration in Fig. 1 has only been used to perform the transformations between a finite group of input/output fields [8,31,32]. In this work, the aim is to calculate the phase masks via the target TM directly. Mathematically, the exact solutions for the phase masks will make the following mean square error (MSE) to zero:

By solving Eq. (4), we obtain a simple analytical solution for φ_{i ,j}. Furthermore, the second-order derivatives of the MSE at these two φ_{i ,j} points can also be easily calculated. It is found that only one of them is positive, which corresponds to the minimum MSE. Hence, there is always only one solution of each phase element φ_{i ,j} for achieving the minimum MSE in one iteration. This is the main finding of our work. The detailed procedure for solving Eq. (4) is presented in Supplement 1. The principle of our algorithm is that the phase elements is calculated one by one. After ramping up all the phase masks once, the program starts from the first phase mask again. The pseudocode of the algorithm is presented in Table 1. After obtaining all the phase values from the algorithm, the phase profiles are discretized to 256 levels for the 0 - 2π values and sent to the SLM.

Table 1. Matrix construction algorithm pseudocode for designing the complex optical transformation matrix realized by cascaded phase masks.

View Table | View all tables in this article

To illustrate the efficiency of our algorithm, we compare our methods with the general-purpose solver for building the optical TM. In Ref. 16, finding a group of phase elements to minimize the MSE, with the constraint of that all phase values lie between −π and π, is treated as a convex problem, which can be expressed as:

In the above equation, the specific properties of the MSE formula are not investigated, thus Eq. (3) is repeatedly calculated in the algorithms. In term of the mathematical operations required, our algorithm exhibits two advantages over the general-purpose solvers.

Firstly, the proposed algorithm can directly obtain each φ_i,j for the minimum MSE with an analytical equation. The general-purpose algorithms are usually based on the framework of gradient descent. The gradient descent computes the gradient of the MSE with respect to each phase element, and changes it along the opposite direction of the gradient. The gradients of the MSEs with respect to every phase element need to be calculated using advanced numerical techniques, as these algorithms are originally designed for the general complicated problems which usually do not have analytical expressions. In contrast, our method studies the concrete optical problem and derives the analytical equation for the first-order derivative of the MSE with respect to each phase element, and confirms that there is only one minimum MSE for each phase value within the range from −π to π.

Secondly, the proposed algorithm avoids abundant matrix operations in Eq. (3) for obtaining the MSEs. During the optimizations, the MSE has to be recalculated once a phase element is updated, so that its gradient can be calculated with respect to another phase element. In the previous algorithm, the matrix multiplications for obtaining the MSEs are the most time-consuming part. With the increases of the number and size of the phase masks, these matrix multiplications become a burden. Whereas, in our method, when solving all the phase values in one phase mask, the MSE can be updated without recalculating the whole TM. The matrix multiplications only need to be performed when one phase mask is completely updated and the optimization for the next phase mask is about to start. To be more specific, when calculating the phase values in the i-th phase mask, we first calculate the TM in Eq. (2) as:

In Eq. (6), T _N,i and T _{i −1,0} are the TMs from the plane after the i-th phase mask to the imaging plane and that from the source plane to the plane before the i-th phase mask, respectively. When one phase value φ_{i ,j} changes, each element in the whole TM (T_pq) can be updated as:

To compare the computing time of the proposed algorithm and the previous methods, we consider a specific optical system with the following parameters: the number of phase mask is 3, the sizes of the phase masks are all 130 × 130, the size of the output field is 16 × 16, and size of the input field is 1 × 1 as being treated as a plane wave. The comparison is performed by using a desktop computer with Intel Xeon Silver 4114 CPU with 20 cores, 256 GB of RAM and the Microsoft Windows 10 operating system. For updating one phase element, the previous general method based on gradient descent derives first the MSE function of the phase element, calculates the differential of the MSE function with respect to that phase element, and finally updates the phase value using the gradient of the MSE and a user-defined learning rate. This process takes about 327 seconds by using our computer. In contrast, our proposed method first calculates the MSE using Eqs. (7) and (S18), and then updates the phase element using Eq. (S17). This process takes about 0.005 seconds under the same computer. The Matlab code for performing the comparison is available in Code 1, Ref. [36], as “method_comparisons.m”.

To the best of our knowledge, this is the first algorithm specially designed for constructing the optical complex-valued TM implemented by the cascaded phase masks. Our algorithm outperforms the existing methods in both computational complexity and efficiency. On one hand, the proposed algorithm reduces the computational complexity by eliminating redundant derivatives and matrix multiplications as aforementioned. The computation complexities of the previous algorithms and our method are compared in Table 2.

Table 2. Computation complexity comparison of the algorithms for obtaining the phase masks via the target complex TM.

View Table | View all tables in this article

On the other hand, the general-purpose algorithms adopt both the gradient of the MSE and a user-defined learning rate to update the phase elements. The learning rate cannot be too large to ensure convergence, while a small learning rate leads to slow convergence. The proposed algorithm calculates the phase value for the minimum MSE in each iteration. Therefore, the decline rate of the MSE in our algorithm is fast. For instance, Ref. 11 updated the matrices over 1000 times in one optimization process, whereas our algorithm reaches to a solution in only 60 iterations, as shown in Table 2.

Apart from the cascaded phase masks, our algorithm has the potential to design other optical architectures, such as the MZI [10] or MMI coupler [11] meshes. However, the following conditions must be satisfied simultaneously: a) the analytical expression exists for the gradient of the MSE with respect to the tuning parameter; b) the MSE function can be effectively calculated and recalculated every time the tuning parameter is updated; c) there is only one value within the tuning parameter’s value range that can achieve the minimum MSE.

Finally, to exclude the possibilities of sticking into local minima of MSEs, we perform multiple runs of the program with random initial phase masks, as it is done in Ref. 11. The results show that the MSEs can always reach to the level of ∼10⁻⁴ regardless of the initial phase values.

3. Numerical results for the construction of the optical linear transform

In this work, we focus on constructing m ² × 1 optical linear complex TM. We numerically implement a common linear transformation, namely, the discrete Fourier transform [11,16] as an example to verify our algorithm numerically. Its expression reads:

By reshaping this 2D function into a vector as described in the previous section, and inserting it into Eq. (3), we obtain the equations to be solved as:

In Eq. (9), we treat the collimated input beam on the first SLM plane as a plane wave with uniform amplitude, and take the TM from the source to the first SLM T ₀ as a vector with all elements equal to 1. For the TM formulated in Eq. (9), the input vector $\vec{X}$ is regarded as a single number to represent the input Gaussian beam. As long as the input beam is collimated, Eq. (9) can be used to simplify the computations without losing generality. The left side of Eq. (9) is the target TM with the dimension of m ² × 1. The right side of Eq. (9) shows the optical TM realized by the multiple reflections from the SLMs. The dimensions of all the phase masks are assumed to be $\sqrt k \times \sqrt k $. As a result, there are in total k elements on each SLM, which can be presented by a diagonal matrix with the size of k × k. The sizes of free-space propagation TM T ₁, T ₂ … T _{N −1} are all k × k, and the size of ${{\boldsymbol T}_N}$ is m ² × k.

In our numerical implementations, the output vector $\vec{Y}$ and each SLM have the size of 256 × 1 (m = 16) and 130 × 130 (k = 130²), respectively. The distances between the neighboring SLMs and that between the last SLM and the imaging plane are all 24.3 cm. In Fig. 2(a) and 2(b), we show the real and imaginary parts of the target TM derived by Eq. (8), respectively. The real and imaginary part of the TM of the constructed optical system consisting of 1 layer of phase mask is shown in Fig. 2(c) and 2(d), respectively. The obtained results are in a good agreement with the target TM. Figure 2(e) compares the MSE converging plots for the optical systems consisting of 1 layer and 3 layers of phase masks. It shows that increasing the number of phase masks can improve the performances of OLT constructions. The 60 iterations for the 1 layer and 3 layers of phase mask(s) take 45 seconds and 1270 seconds by using our computer, respectively. More detailed comparison of the reconstructed TM realized by 1, 2 and 3 layer(s) of phase masks is shown in Fig. S1 in Supplement 1.

 

Fig. 2. The numerical results of the optical complex TM designed by the proposed algorithm. (a) and (b) The real and imaginary parts of the target optical complex TM. (c) and (d) The real and imaginary parts of the numerically implemented optical complex TM constructed by the optical system consisting of 1 layer of phase mask with the size of 130 × 130. (e) The convergence plots of the MSEs over the iterations for the optical systems consisting of 1 layer and 3 layers of phase masks with the size of 130 × 130. Here, one iteration includes the updates of all phase elements once.

Download Full Size | PDF

To further illustrate the performances of the proposed algorithm for optical configurations with different numbers of phase masks, we compare the achieved MSEs and computing time for different layers of phase masks in Table 3. From the numerical results, it can be concluded that the optical configurations with more phase masks can achieve lower MSE if the numbers of phase elements in one phase mask and the iteration numbers are kept the same. However, increasing the number of phase masks will also increase the computing time. It should also be noted that the reductions of MSEs for the optical configurations with more phase masks are larger after the same number of iterations. So, the optical configurations with more phase masks actually require less iterations. For instance, it takes only 34 iterations for the 3 layers of phase masks to achieve the same MSE that the single layer phase mask achieves after 60 iterations, as shown in Fig. 2(e). The Matlab code for performing these numerical simulations is available in Code 1, Ref. [36], as “phase_retrieval_for_DFT_in_simulations.m”.

Table 3. Comparison of simulation results for different numbers of phase masks

View Table | View all tables in this article

To choose the number of phase masks in the configurations, one must make the trade-off between the benefit of achieving low MSEs and the drawbacks of increasing the optimization time and the complexity of the experimental setup. The appropriate number of phase masks should be determined by the sizes of phase masks, input and output optical fields. For the applications considered in this work, we found that 1 layer of phase mask is sufficient for delivering low-MSE results and can lead to a much simpler experimental setup.

4. Experiment results for the generation of optical complex data

As concluded in the previous section, we use only one layer of phase mask in our experiments in order to keep the whole optical system simple, as depicted in Fig. 3. The detailed experimental setup is explained in Fig. S2 in Supplement 1. The optical system in Fig. 3 can be formulated as the right side of Eq. (9) by taking N = 1. By setting a target TM as the left side of Eq. (9) and employing the proposed algorithm to solve Eq. (9), an optical array can be experimentally generated. Each beam is an optical manifestation of a complex number. In this section, we experimentally demonstrate that the beams in the optical array are independent, and their amplitudes and phases can be separately controlled.

 

Fig. 3. The schematic diagram of the experimental set-up used to verify the proposed matrix construction algorithm. The experiment generates the complex-valued optical array by phase-only modulation. The optical system is composed of a free-space propagation from the source to the SLM (T ₀), an SLM (diagonal matrix containing phase delay e^iφ) and a free-space propagation from the SLM to the imaging plane (T ₁). The optical system transforms a solo input beam into an array of 64 beams and each beam corresponds to an independent complex number. The distance between the SLM and the CCD planes is 33.5 cm.

Download Full Size | PDF

Before we implement the proposed method to generate optical complex data, we first need to address the uniformity issue, a common problem in the multi-beam generation scenarios, as an essential preparation step. Due to the non-uniform intensity distribution of the incident beam and the imperfect phase implementations, the intensities of the beams away from the array center are smaller than those in the central part, even though they are designed to be the same. This problem also causes distortion and loss of information in the computer-generated hologram (CGH) techniques. This imperfection can be compensated by iteratively setting a new target intensity distribution according to the resulted images in the previous iterations when generating the CGHs [37–39]. Here, we show that our algorithm can also compensate for the observed intensity non-uniformity. The intensity calibration is implemented by adjusting the intensities in the target TM according to the actual experiments, i.e., the designed intensities of the experimentally obtained strong and weak beams are decreased and increased in the next matrix construction process, respectively. The CCD images of the 8 × 8 beam arrays before and after the intensity calibration are shown in Fig. 4(a) and 4(b), respectively. Figure 4(c) shows the digital data from the CCD images, summed along the vertical direction within the region denoted by the red dashed line. The results show that our algorithm can yield uniform beam arrays. The calibration results of the arrays composed of 4 × 4 and 8 × 8 beams and the corresponding phase masks are shown in Fig. S3 and S4 in Supplement 1, respectively.

 

Fig. 4. The intensity calibration results of the beam arrays. The CCD images of the 8 × 8 beam arrays captured (a) before and (b) after the intensity calibration. (c) The digital data of the CCD images, summed along the vertical direction within the region denoted by the red dashed lines, plotted along the red dashed lines in (a) and (b).

Download Full Size | PDF

To demonstrate that the intensity of each beam can be independently controlled through our method, we define multiple target TMs with different intensities of one parameter and feed them into the proposed algorithm to retrieve the phase masks. Figure 5 shows the experimental results for varying the designed intensity of one chosen beam in the target TM, indicated by the red square. The results in Fig. 5(a) are obtained by averaging the digital data on the CCD image over the red square region. Ideally, the relation between the intensity of the generated beam and the corresponding designed intensity should be linear. However, when the intensity of one beam is low, the CCD cannot well display the image of that beam. Therefore, in the practical experiment, the measured result is nearly linearly proportional to the designed intensity only when it is larger than 0.2. This relation in Fig. 5(a) can be fitted as:

Figure 5(b) and 5(c) show the CCD images when the amplitude of the beam indicated by the red square is set to 0.6 and 1, respectively. It clear that the marked beam intensity can be manipulated, while the rest beams are barely affected. The results in Fig. 5 are obtained after the aforementioned intensity calibration. More details on varying the intensities of the selected beams are shown in Visualization 1. Naturally multiple beams can also be controlled simultaneously, as other TMs can also be set as the targets in the proposed algorithm.

 

Fig. 5. The experimental results for manipulating the intensity of one beam in the optical array. (a) The relation between the average digital data captured by the CCD over the red square region and the designed intensity. (b) The CCD image showing the beam arrays with the designed intensity of the beam in the red square set to 0.6, while the others set to 1. (c) The CCD image showing the beam array with the designed intensities of all beams set to 1.

Download Full Size | PDF

To illustrate the phase control of beams, we use a second SLM as a checker (see Fig. S2 in Supplement 1). The checker SLM steers one beam within the array and makes it interfere with other beams, which leads to fringe patterns. The phase profiles of the three beams passing the checker SLM and the intensity distributions of the three beams without the checker SLM are shown in the upper and lower row of Fig. 6(a), respectively. As a result, by observing the shift of the fringes, we can identify the phase difference between two beams. The fringe patterns from the three beams with phases of θ ₁ = 0°, θ ₂ = 0°, θ ₃ = 36° and θ ₁ = 0°, θ ₂ = 0°, θ ₃ = 216° are shown in Fig. 6(b) and 6(d), respectively. As θ ₃ in the target TM is changed, the corresponding fringe shifts, meanwhile, the other interference pattern is barely affected. To identify the change of θ ₃ through the fringes, the profiles of the CCD results along the horizontal direction is plotted in Fig. 6(c) and 6(e). The CCD results are the summation of the digital data from the CCD images in Fig. 6(b) and 6(d) along the vertical direction. The difference between the two values θ ₃ in Fig. 6(b) and 6(d) can be calculated as:

In Eq. (11), x ₂₂−x ₁₁ is the period of the phase profile of the steered beam from the middle beam. This result is consistent with the design. It is worth noting that the pixel size of the CCD is 0.02 mm, and the position error of even one pixel will cause significant deviation for the calculated phases. A more precise phase measurement could be conducted with a Shack–Hartmann wavefront sensor [22]. The different shapes of the two fringes may result from the uneven amplitude splitting of the middle beam. The different intensities of the middle beams in Fig. 6(b) and 6(d) may be caused by the potential partial overlap of beams on the checker SLM. Despite the distortion and noises, the results indicate that our algorithm can indeed generate independent phase values among the beams. More details on varying the phase of one of the beams are shown in Visualization 2. By resetting multiple elements in the target TMs, the phases of multiple beams can also be controlled simultaneously, hence various complex-valued vectors can be encoded into the optical arrays. The Matlab code for generating the phase masks used to obtain the experimental results is available in Code 1, Ref. [36], as “phase_retrieval_for_beam_array_in_experiments.m”.

 

Fig. 6. The experimental results for manipulating the phase of one beam in the optical array. (a) The three beams with independent phases: the phase profiles of the three beams passing the checker SLM and the intensity distributions of the three beams without the checker SLM are shown in the upper and lower row, respectively. The beams passing through the checker SLM can produce two interference fringes. (b) and (d) The fringe patterns from the three beams with phases of θ ₁ = 0°, θ ₂ = 0°, θ ₃ = 36° and θ ₁ = 0°, θ ₂ = 0°, θ ₃ = 216°, respectively. (c) and (e) The digital data of the CCD images, summed along the vertical direction, plotted along the horizontal direction in (b) and (d), respectively.

Download Full Size | PDF

Finally, we compare our scheme with other optical beam array generation technologies, e.g., phase grating and CGH. The phase grating method needs a Fourier lens to image the beam arrays at its focus plane. This scheme can also control the amplitude and phase of each beam, but changing the amplitude of one beam will to a large extend affect other beams as well. Another disadvantage of the phase grating method is that the undesired high diffraction orders are usually inevitable [24]. In contrast, our method can obtain the aforementioned results without using Fourier lenses nor spatial filters, and the generated optical arrays show negligible unwanted beams beyond the designed arrays. On the other hand, the CGH methods use various iterative algorithms, e.g., Gerchberg–Saxton (GS) algorithm [39], to retrieve phase mask for a target beam array. However, the phase of each beam cannot be controlled, as the phase information is a free parameter in the GS algorithm. In comparison, our method can manipulate both amplitude and phase of beams individually. The comparison of experimentally realized beam arrays using the above methods is shown in Fig. 7.

 

Fig. 7. Beam arrays experimentally realized with different methods. (a) phase grating, (b) the GS algorithm for generating CGH, (c) the proposed matrix construction algorithm (this work). The red box in (a) indicate the beams that are intended to be generated.

Download Full Size | PDF

Moreover, the numerical spatial frequency filter, which has been previously applied in some CGH algorithms [23,34], can also be incorporated into our algorithm to eliminate the interferences from the high spatial frequency components, resulting in a smooth image without adding any optical elements. The detail implementation of the numerical spatial frequency filtering and the corresponding results are described in Section 5 in Supplement 1.

5. Conclusion

In this work, we have proposed the first algorithm specially designed for constructing the optical linear complex matrix. As a direct application of our algorithm, we have implemented an optical complex data generator using phase-only modulations. Firstly, we experimentally demonstrate that the proposed algorithm can compensate the imperfections in the light path, thereby creating beams with uniform intensities. This is an essential preparation step for encoding multiple complex numbers into these beams. Next, the experiments confirm that every beam in the optical arrays is independent, and their amplitudes and phases can be individually controlled. Compared with other optical beam array generation techniques, our scheme can suppress the unwanted noises, offer the full flexibility to control both the amplitude and phase of each individual beam, and be readily incorporated into most current optical computing architectures thanks to the simple optical set-up consisting only of a light source and one phase mask.

Apart from the aforementioned application, the significance of this work is the discovery that there is an analytical solution for each phase value in the OLT system, which is derived from the equations connecting the formula for the light propagation through phase masks and the given complex-valued matrix. This direct calculation method is much more efficient than the commonly adopted deep learning-based methods [9]. Our algorithm also exhibits much higher efficiency than the general-purpose algorithms, as the redundant matrix multiplications and differential calculi in every iteration can be eliminated altogether. Our work provides a new insight into the design of the general OLT systems.

Although only the transformation between a Gaussian beam and an optical beam array described by the m ² × 1 complex matrix has been experimentally demonstrated here, the proposed algorithm can also be used to realize the more general m ² × n ² matrix. It must be pointed out that our algorithm addresses the design of such an optical system through its TM directly, not from a finite group of input/output fields (the n ² × 1 and m ² × 1 vectors). Therefore, besides its high efficiency, another important advantage of the proposed algorithm is that the constructed OLT is irrelevant to the input light fields and can perform the same all-optical linear computations to any given n ² × 1 input vector. In contrast, the solutions derived from a set of input/output fields only approximate the required TMs, and the resulted all-optical operations are less accurate when applied to the input fields outside the design-fields group [40].

The proposed optical complex data generation method still suffers from two major problems for practical applications. Firstly, the data processing speed of the optical computing scheme based on the commercial phase-only SLMs is still limited by the refresh rate of liquid crystals. It is expected that the optical complex data can be generated or computed at much higher speed by employing the recently demonstrated fast photothermal SLMs [41]. Secondly, the overall size of the optical setup is still bulky, as the free-space optical propagation is critical to allow light to diffract and interfere to construct a desired optical pattern. The recently proposed “spaceplate” concept has the potential to shrink the future optical computing systems by squeezing free space [42].

Our future work will be dedicated to the deeper exploitation of the proposed algorithm for more advanced all-optical functionalities. For example, by further increasing the constructed matrix dimensions, our algorithm may find broader applications in building large-scale optical parallel processing networks. Currently, the proposed design strategy is only applicable to the optical architecture with the cascaded phase masks. We will analyze the mathematical expressions for the other OLT implementing techniques summarized in [43] to find whether our design strategy is also applicable, or develop specific design procedures in our future work.