1. INTRODUCTION

The ability to perform linear operations on light with optical devices is a fundamental ingredient in many areas of optics and photonics, including signal processing and spatial multiplexing in optical communication, as well as optical analog computation. For emerging applications to optical artificial neural networks, the ability to reconfigure these linear photonic processing units is crucial. Moreover, ease of fabrication and robustness to manufacturing inaccuracies are highly sought-after features. Similar to electronic devices, photonic devices are traditionally designed to perform one given operation [1–5]. The conformation of the device is directly linked to its intended function. As a consequence, fabrication imperfections and changes of environmental conditions negatively affect its functioning, limiting the range of operation of the device. Furthermore, such inverse-design approaches inherently prohibit reconfigurability.

Programmable coherent photonic circuits have the potential to offer reconfigurability. Initially conceived in free-space optics by leveraging beamsplitters [6], recent advances in silicon photonics also enabled the implementation of the concept in integrated designs [7–12]. Several first experimental demonstrations of this approach have produced promising results. However, the amount of phase shifters required in these architectures scales with the square of the size of the desired transformation, which raises questions about the overall scalability of the approach. Alternatively, the idea of controlling mode coupling inside a multimode waveguide to shape the output wavefront was demonstrated using silicon-on-insulator planar waveguides whose index profile can be tuned using femtosecond laser pulses [13]. Another approach using multiplane modulation was proposed [14] and was successfully used for few mode manipulations [15]. While small transformations may suffice for optical communications, the increasing complexity of the physical system may hinder the implementation of the high rank transformations needed in optical analog computation. There is currently no solution to perform optical linear operations involving a high number of channels.

Parallel to the aforementioned developments, since Vellekoop and Mosk’s landmark paper in 2007 [16] the field of wavefront shaping in complex media emerged [17,18]. A complex medium may be defined as a system that mixes the spatial and/or temporal degrees of freedom, resulting in the complete scrambling of an incident wavefront [19]. Examples include chaotic cavities, disordered waveguides, or random scattering systems such as paint layers or biological tissues [16,17,20–25]. The seemingly random effect of the disorder on the wavefront is deterministic; hence, such a system can be fully represented by a linear transmission matrix (TM) [20]. Due to the large number of modes, those matrices have large dimensions. Moreover, the lack of symmetry of the system contributes to the TM’s high rank. By injecting light into a well-chosen combination of input modes using a spatial light modulator (SLM), an optical scattering medium can be used to perform a wide variety of functions. Initial efforts sought to overcome the scrambling of the wavefront in space and/or time through wavefront shaping (e.g., to focus behind or transmit information through a complex medium) [16,20,25–30]. Strikingly, it was possible to beat the Rayleigh limit on focusing [31,32] and the Nyquist–Shannon sampling theorem [33].

Despite its potential, controlling wave propagation in complex media to perform linear transformations has received little attention. The literature contains reports on implementations of two-port beamsplitters [34,35], with applications to control quantum interferences in mind [36,37], as well as on spatial mode sorting [38]. In particular, the possibility to implement large reconfigurable linear transformations as required for optical neural networks has remained unexplored to date. Building on recent work that demonstrated wave-based analog computation in a chaotic microwave cavity [39], we explore here the possibility to perform large complex-valued linear operations in optical complex media. We formulate fundamental upper bounds in terms of the operation size that can be implemented in a realistic optical system. Moreover, we experimentally demonstrate the possibility to take advantage of the large number of degrees of freedom of a disordered optical system (multimode fiber or glass diffuser) of which we acquire the transmission matrix (TM), to physically perform different desired linear transformations. We aim to provide an implementation of reconfigurable optical linear transformations that meets the necessary criteria for real-life implementations (i.e., minimal hardware requirements, low stability requirements, and short reconfiguration times). To that end, we do not rely on interferometric phase measurements; instead, we leverage phase retrieval techniques, eliminating the need for high mechanical stability. We use fast binary modulators that allow above 10 kHz refresh rates. Moreover, the one-off calibration of the medium’s TM allows us to reconfigure the implemented linear transformation without any further experimental measurement.

2. OPERATION PRINCIPLE

A. Principle

Light propagation through a linear disordered system at a given frequency is fully described by its transmission operator $\mathbf{H}$ that links the output state of light $|{\psi }_{\mathrm{out}}⟩$ to the input one $|{\psi }_{\mathrm{in}}⟩$, so

B. Input and Output Projectors

We limit ourselves to input projectors ${\mathbf{P}}_{\mathbf{in}}$ that can be created by modulating the optical field at a given input plane (i.e., using a SLM). We divide the SLM into $n$ groups of $N/n$ pixels on which we control the amplitude and/or the phase of the optical field. The output projection ${\mathbf{P}}_{\mathbf{out}}$ is realized by measuring $m$ speckle grains at the output plane.

For illustration purposes, we use $n=m=4$. The corresponding projectors ${\mathbf{P}}_{\mathbf{in}}$ and ${\mathbf{P}}_{\mathbf{out}}$ have the following matrix representations,

C. Existence of Solution with Ideal Modulation and Theoretical Limits

We assume the ideal case of $\mathbf{H}$ being a Gaussian random matrix measured with negligible noise and the spatial modulation scheme being able to control with high fidelity the amplitude and the phase of the field. The Gaussian matrix approximation is known to be valid when one controls and measures a fraction of the total number of modes [40]. Satisfying the equality of Eq. (2) with the projectors as represented by Eq. (3) consists in solving $n$ systems of $m$ linear equations, so

Note that no hypothesis has been made on the operator $\mathbf{G}$. It demonstrates the universality of the approach to perform any linear transformation. Experimentally, the passive modulation of the SLM imposes $|p_{i,j}|<1\forall (i,j)\in [1,N/n]\times [1,n]$. It limits the applications to operators with below unity singular values. However, for any operator $\mathbf{G}$, the experimental system we present can simulate the transformation corresponding to the operator ${\mathbf{G}}^{\prime }=\alpha \mathbf{G}$ with $\alpha \le 1$ chosen so that the largest singular value of ${\mathbf{G}}^{\prime }$ is smaller than 1.

3. EXPERIMENT

A. Optical Setup

To illustrate the versatility of our approach, we report two experimental implementations of our scheme using two very different complex media as optical processing units: a multimode fiber and a ground glass diffuser at two different wavelengths (respectively, 1550 nm and 632.8 nm). For both experiments, we modulate the impinging wavefront with a digital micromirror device (DMD).

The MMF experiment uses a 1 m long segment of a step index fiber of numerical aperture (NA) 0.22 with a core diameter of 105 μm supporting approximately 2000 guided modes. The light source is a 1.55 μm telecom narrow-band laser (PS-NLL, TeraXion). We use a $1920\times 1200$ pixel, 16 kHz digital micromirror device (V-9601, ViALUX GmbH) for the modulation. The output intensity pattern is imaged onto a fast $640\times 512$ pixel InGaAs camera (Cheetah 640CL, Xenics). The fiber is compressed at four different locations to ensure strong mode coupling [43]. The scattering media setup uses a ground glass diffuser (DG20-1500, Thorlabs) with a 632.8 nm laser source (1137/P, JDSU Uniphase) modulated by a $2560\times 1600$ pixel, 12 kHz DMD (V-9001, ViALUX GmbH) and a $2336\times 1752$ pixel CCD camera (Prosilica GT2300, Allied Vision). The MMF setup is depicted in Fig. 1, and the scattering medium setup is schematically similar. The number of controlled segments on the DMD is $N=1568$ for the MMF experiment and $N=2304$ for the scattering media. Both complex media stay highly correlated during the time of the experiment: The output intensity patterns show correlations greater than 98% for about 2 h. (Supplement 1 provides the details and correlation curves.)

 

Fig. 1. Schematic representation of the setup. An expended laser beam is modulated after reflection off a DMD and injected into a complex medium (CM: ground glass diffuser or multimode fiber). One polarization of the outgoing light is recorded by a digital camera. A reference arm is used only for the final estimation of the fidelity to measure the complex optical field for the MMF experiment. FC: fiber coupler; Ci (with $i\in [1\dots 4]$): fiber collimator; Li (with $i\in [1\dots 4]$): planoconvex lens; PBS: polarization beamsplitter; and P1 and P2: polarizers.

Download Full Size | PDF

B. Procedure

The steps of the procedure are schematically illustrated in Fig. 2. A first calibration step consists in learning the TM of the complex propagation medium, namely a multimode fiber or a scattering medium. The second step consists in numerically finding the optimal projectors ${\mathbf{P}}_{\mathbf{in}}$ and ${\mathbf{P}}_{\mathbf{out}}$ that best satisfy Eq. (2) for a given operator $\mathbf{G}$ knowing the transmission matrix $\mathbf{H}$. Finally, by displaying the corresponding masks onto the SLM and measuring the field on the designated outputs, the system acts as the target operator $\mathbf{G}$ on the incident field. Those three steps are discussed in detail in the following sections.

 

Fig. 2. Overview of the experimental procedure for performing a DFT operation of size 4 ($n=m=4$ and $\mathbf{G}={\mathrm{DFT}}_4$). Step 1: System calibration. Acquisition of the complex TM by measuring a set of input patterns and output intensity speckles and using a phase retrieval algorithm. Step 2: Optimal input projection calculation. We illustrate the procedure by showing how to find one subpart of the input mask. The computation is done independently for each subpart of the SLM using a convex optimization solver and using the TM and the corresponding column of the target matrix $\mathbf{G}$ as inputs. Step 3: Analog computation. The proposed optical processing unit is composed of a spatial light modulator (SLM) and a complex medium. The SLM and the output detection take the role of the projectors ${\mathbf{P}}_{\mathbf{in}}$ and ${\mathbf{P}}_{\mathbf{out}}$, converting the given transmission matrix $\mathbf{H}$ of the complex medium into a desired linear transformation $\mathbf{G}$ [see Eq. (2)].

Download Full Size | PDF

1. Calibration

We first estimate the medium’s complex-valued transmission matrix $\mathbf{H}$. To remove the need for interferometric measurements, which requires high mechanical stability that limits the versatility of the approach, and to mitigate the effect of measurement errors, we use a phase retrieval algorithm. It allows the recovery of the full TM from intensity-only measurements [44,45]. We divide the pixel array of the DMD into $N$ groups of pixels or macropixels. We then send a learning set of $7N$ random binary vectors (entries have zero or unity amplitude) and measure the corresponding output intensity patterns on the camera. These patterns are projected onto an output basis of macropixels of a size that corresponds to a speckle grain. The obtained data is fed into a phase retrieval algorithm [45] and the results are further refined using a gradient descent optimization to better account for the nonlinearity of the camera response (see Supplement 1).

We reconstruct matrices of size $100\times 1568$ for the MMF setup and $100\times 2304$ for the scattering medium. The computation is accelerated through parallelization with a graphical processing unit (GPU) (GeForce GTX 1050 Ti, NVIDIA) and is completed in under 2 min. We evaluate the quality of the reconstructed TM using test input wavefronts that have not been used for the learning procedure by calculating the root mean square error (RMSE) for each element of the output basis, giving ${\mathrm{RMSE}}_{\mathrm{mean}}\sim 7\%$ ($\mathrm{std}({\mathrm{RMSE}}_{\mathrm{mean}})\sim 6\%$ and ${\mathrm{RMSE}}_{\text{median}}<5\%$) for the MMF setup, and ${\mathrm{RMSE}}_{\mathrm{mean}}\sim 11.6\%$ ($\mathrm{std}({\mathrm{RMSE}}_{\mathrm{mean}})\sim 7.3\%$ and ${\mathrm{RMSE}}_{\text{median}}<7.9\%$) for the visible setup.

2. Projection Calculation

The output projection ${\mathbf{P}}_{\mathbf{out}}$ is realized by selecting $m$ output speckle grains. The optical fields in these areas correspond to the output of our operator. We choose the output points that give the lowest reconstruction error during the calibration step. We then obtain the output projection as represented in Eq. (3), where rows are ordered from the lowest to the highest error.

We demonstrated that for a full amplitude and phase modulation, the input masks corresponding to the optimal input projectors ${\mathbf{P}}_{\mathbf{in}}$ can readily be calculated from the TM using the relation 5. Using binary amplitude modulators, we achieve a few-level modulation of the optical phase using Lee holograms (see Supplement 1). The modulated beam is then projected onto the complex medium. We collect the output intensity pattern for a single polarization using a digital camera. Identifying an input mask on the DMD that approximates the equality in Eq. (2) for a given target matrix $\mathbf{G}$ is an ill-posed problem since we do not have full control over the complex wavefront. Finding an input projector ${\mathbf{P}}_{\mathbf{in}}$ that minimizes $\Vert \mathbf{G}-{\mathbf{P}}_{\mathbf{out}}^{\mathbf{T}}{\mathbf{HP}}_{\mathbf{in}}{\Vert }_2$ with the constraint of looking for solutions with discrete phase values is a convex problem. We use a mixed-integer convex solver [46,47] to find an approximate solution that satisfies the experimental restrictions of the few-level wavefront modulation of the DMD (see Supplement 1). The computation of the input projector ${\mathbf{P}}_{\mathbf{in}}$ takes less than 20 min for $16\times 16$ operators and under 5 min for $8\times 8$ operators on a computer with an Intel i7-7700 CPU with 32 GB or RAM running on Windows 10.

3. Optical Analog Computation

Once the appropriate input projector mask is calculated and displayed on the DMD, the system is ready to act like the desired linear operator $\mathbf{G}$. The operation is performed on a desired input vector by encoding its components $\mathbf{X}={(x_1,x_2,\dots ,x_n)}^T$ into as many light beams directed toward the SLM and collecting the output vector $\mathbf{Y}$ on the detection device such that

4. RESULTS

We illustrate the reconfigurability of the presented scheme by experimentally implementing two linear transformations very common in computer and physical science; namely, the discrete Fourier transform and the Hadamard matrix. Their general expressions read

To measure the complex output field, we carry out off-axis holography measurements [48]. It uses a reference arm originating from a 90/10 fiber splitter at the output of the laser. The reference beam and the output of the MMF are recombined using a polarizing beam splitter and a polarizer. The recorded images are post-processed to retrieve the phase information from the interference pattern. Note that the reference arm was not part of the procedure of implementing the linear operator; it is solely used to monitor the effective complex-valued operator $\tilde{\mathbf{G}}$ afterward. (Supplement 1 provides details about the complex field measurements.) Moreover, as a single image allows measuring the complex value output of the operator, interferometric stability over times longer than the integration time of the camera is not required.

We experimentally perform all optical operations according to the presented principle. A large range of input signals is prepared and sent to the setup (see Supplement 1). The quality of our results is assessed by the Pearson correlation coefficient between the absolute values of the experimental and the ideal output vectors $\mathbf{C}$ and also by the experimental fidelity ${\mathbf{F}}_{\mathbf{c}}$.

In Figs. 3(a) and 3(c), we compare the measured amplitude values obtained for two random input vectors with the expected ones for both operators. The results show good agreement, as emphasized by the correlation coefficients exhibiting values over 95% for these realizations. In Figs. 3(b) and 3(d), we show the real and imaginary parts of the estimated experimental operators $\tilde{\mathbf{G}}$ after normalization and the real and imaginary parts of the target operator $\mathbf{G}$. The results presented were obtained without averaging. A summary of the measured quality estimators for the MMF experiment is presented in Table 1 with and without averaging for operator sizes of $n=m=8$ and $n=m=16$.

Table 1. Summary of the Efficiency Results for the MMF Experiment^a

View Table

 

Fig. 3. Comparison of experimentally implemented operators and target operators. We use as target matrix $\mathbf{G}$, a $16\times 16$ Hadamard matrix [(a) and (b)], and a $16\times 16$ discrete Fourier transform [(c) and (d)]. Absolute values of experimental (red dots) and theoretical (cyan diamonds) results of the linear transformations for two different random input vectors drawn from $\{-1,0,1\}$ for $\mathbf{G}={\mathrm{Ha}}_{16}$ (a) and $\mathbf{G}={\mathrm{DFT}}_{16}$ (c). The correlation between the two signals is shown in the insert. Comparison between experimental and target operators for $\mathbf{G}={\mathrm{Ha}}_{16}$, $F_c=0.785$ (b) and $\mathbf{G}={\mathrm{DFT}}_{16}$, $F_c=0.818$ (d). The results are obtained without averaging.

Download Full Size | PDF

A good agreement between the experimental data and the ideal operations outputs is observed for the different sizes tested, even without averaging. It demonstrates the possibility of performing one shot operations through a multimode fiber. To further illustrate the possibility to create any desired operator, in Fig. S5 (Supplement 1) we show an experimental implementation of a matrix with no physical meaning displaying the name of the authors’ host institution by encoding information independently in the real and imaginary part of its elements.

Because the ground glass diffuser experiment is more prone to errors in the operator reconstruction, we only present results for $n=m=8$. In Fig. 4 we show the raw amplitude of four different outputs for the ${\mathrm{DFT}}_8$ operator. We measure $\mathbf{C}=0.948\pm 0.018$. We estimate that the sources of the errors are the lower laser stability and the nonlinear response of the CCD camera. While negatively affecting the results, we obtain qualitatively similar results, demonstrating the capacity of the approach to work with relatively cheap devices, namely a HeNe laser and a standard CCD camera. Results further illustrating the versatility of our method are displayed in Supplement 1.

 

Fig. 4. Amplitude of four different random output vectors for operator $\mathbf{G}={\mathrm{DFT}}_8$ obtained using a ground glass diffuser as a complex medium. Cyan diamonds correspond to the theoretical outputs and red dots to the experimental data. Inserts in the top left corners give the correlation between the experimental and the predicted transforms.

Download Full Size | PDF

5. DISCUSSION

We demonstrated the possibility of using cheap and common media as optical processing units to perform linear operations using wavefront shaping. The attractiveness of such an approach is linked to its ability to be scaled up for larger operations. While we report the implementations of $16\times 16$ operators, much larger operations can be performed, provided there is an increased control over the input field and a reduced noise level. Noticeably, Eq. (2) was only satisfied approximately because of the binary nature of the DMD modulation; the search for such approximation requires computation efforts. While averaging over realizations mitigates the error, a single shot operation is usually desirable. For an ideal amplitude and phase modulation, Eq. (2) can be satisfied exactly with a simple matrix pseudo-inversion. A complete independent amplitude and phase modulation can for example be obtained using two liquid crystal SLMs. Another approach is, after characterization of the medium, to print phase plates to fix the input mask. While limiting the reconfigurability of the system, it reduces the final cost of the processing unit.

For a square target matrix $\mathbf{G}$, the largest size for which one can find a solution for a fixed number $M$ of controlled modes is $n=m=\sqrt{M}$. Using a multimode fiber as a complex medium, the limiting factor is the number of modes supported by the fiber; for a typical large core step index fiber (550 μm core, $\mathrm{NA}=0.22$), $M\approx 60,000$, corresponding to maximal operator dimensions $n=m\approx 250$. In contrast, in scattering layers such as glass diffusers, the number of degrees of freedom available given by the number of propagating modes is quasi unlimited. Thus, the number of independently controlled modes is limited by the number of pixels on the modulator, typically on the order of one million. Note that this remains true in the multiple scattering regime as long as the number of transmission channels is large compared to the number of input pixels (i.e., far away from the localization regime). Hence, this would allow creating linear operators of size $n=m\approx 1000$. These operator sizes match the order of magnitude of the size needed for optical neural networks.

It is important to note that our approach, due to intrinsic losses and the absence of gain in the system, can perform any linear operation only up to a constant multiplier. In particular, it cannot perform operations with above unity singular values. Such a restriction can be detrimental to quantum optics applications where losses can modify the optical state of light. Our apparatus based on DMDs causes more than 50% of the light to be lost upon modulation. However, using a phase-only SLM based on deformable mirrors or a phase plate together with a careful injection into a multimode fiber, close to unitary operations can be achieved. Multimode fibers have the potential to outperform integrated-photonics-based platforms in terms of losses, with insertion losses on the order of 0.3 dB and propagation losses below 1 dB/km, compared to values above 1 dB and 0.1 dB/cm, respectively, for photonic integrated circuits [49].

The presented method requires computational efforts during its calibration; however, once the TM is retrieved and the projectors are calculated, the operations are performed in a single shot ($O(1)$ operations) on a passive system. No further calibration is needed as long as the system is stable. Multimode fibers up to 100 m are stable enough to be used as accurate optical instruments [50]. Stability over multiple days can be obtained by ensuring a controlled environment and good laser stability. Thin diffusers are only limited by the pointing stability. The proposal offers the possibility to implement large optical linear transformations without elaborate fabrication techniques as well as to reconfigure the desired operator without further measurements. Moreover, it opens the opportunity to drastically reduce the energy consumption compared to classical electronic components while increasing the computation speed [11,39]. These characteristics may enable the presented technique to play a key role in the advent of optical analog computation and machine learning.