Emulating quantum computing with optical matrix multiplication

Controlling various degrees of freedom of light, i.e., time, frequency, space, and momentum, has become an emerging and promising tool for numerous information processing tasks in classical and quantum domains, ranging from novel imaging methods^1–7 to communications^8–13 to computations,^14–17 all harnessing the high dimensional nature of structured light fields.^18–22 This is because spatial modes enable information encoding in qudit spaces for dimensions d > 2 per particle as opposed to the d = 2 (qubits) encoding levels offered by traditional qubit encoding, such as with polarization states. Because qubits are easy to control, they are a reliable resource for numerous protocols.

In traditional quantum computing, encoding basis states are typically formed from multiple qubits, and the operations needed to construct arbitrary unitary gates in d > 2 dimensions often involve a mix of two-dimensional gate operations (such as the Hadamard gate and T-gate) along with multi-qubit operations (such as the Control-not gate).²³ This allows for precise and efficient execution of arbitrary computations that are required to demonstrate quantum advantage.²⁴ Given that each qubit is restricted to two dimensions, the only way to enhance encoding capacity is by increasing the number of qubits and ensuring optimal connectivity among them. An alternative is to utilize qudits, which are particles that occupy d > 2 dimensions because they offer larger encoding state spaces²⁵ and promise robustness against noise.²⁶ As a result, structured photon states emerge as a promising candidate for achieving this, as they occupy higher-dimensional Hilbert spaces. Furthermore, it has been demonstrated that vector-matrix multiplication, a fundamental operation in quantum mechanics, can be performed using transverse optical fields.²⁷ 

There have been numerous impressive proposals and implementations of quantum computing algorithms that utilize the higher dimensional nature of the internal degrees of freedom of photons.^28–30 This field has developed tremendously, aided by the advent of multiplane light conversion technology,³¹ where digital holography and free-space diffraction are employed to realize arbitrary unitary gate operations.³² Inspired by these advances, diffraction-based deep neural network architectures have also been proposed.³³ In addition, harnessing mode mixing in complex media such as multimodal optical fibers, as opposed to diffractive MPLCs, has emerged as a potential avenue (see review on this topic³⁴) to achieve similar operations but with discrete pixel-like states.³⁵ To realize the full potential of qudits for quantum computing, efficient single photon sources and the tools for entangling and controlling them are necessary.

On the other hand, harnessing classical optical fields for quantum computing has also been an active area of research.^14,36–40 The drive in this field has been inspired by the observation that some of the essential resources of quantum computing can be found in classical coherent fields. Similar to quantum states, classical waves can be prepared in superpositions and can be interfered with, allowing for parallel information processing. Notably, Jozsa also argued that classical waves can provide an efficient simulation platform since not all quantum algorithms require entanglement as a resource.^41,42 Recent advances have demonstrated that optical methods can implement unitary and nonunitary transformations efficiently, utilizing programmable holographic techniques and spatial phase modulation to enable high-dimensional operations.⁴³ Similarly, parallel processing for computational tasks such as multiplication modulo, a crucial component of Shor’s algorithm, using phase modulation has been explored, offering promising avenues for structured light applications in computational tasks.⁴⁴ This opens up new possibilities for studying or implementing quantum computational tasks without relying on complex quantum hardware by merely exploiting coherent fields and performing the required matrix-vector operations on them.

An intriguing and unexplored approach to utilizing coherent fields in quantum computing is Tamura’s optical matrix-vector multiplication from the 1970s,⁴⁵ which was long forgotten but has recently resurfaced. This method harnesses point-wise multiplication (Hadamard product) of light and the Fourier transform capabilities of lenses. With this approach, operations have been applied in dimensions reaching up to d = 56 dimensions.²⁷ Hidden in this method is the capability to encode information in Hilbert spaces formed from transverse modes of coherent laser fields.

In this work, we use optical matrix multiplication as a foundation for emulating quantum algorithms, therefore interacting quantum computing with structured light. We unveil the inherent tensor product structure of transverse spatial modes that is exploited in coherent optical matrix-vector multiplication and show that it can be harnessed as a tool to emulate quantum computing algorithms. Our logical basis is constructed from lattices of displaced Gaussian modes, eigenmodes of free space, that are locally modulated with digital holograms encoded on spatial light modulators. We treat the Cartesian location coordinates x and y of each mode as our degrees of freedom, forming the equivalent of two qudit registers. The x component is used to embed our state vectors, while the y-component acts as an ancilla that facilitates the unitary gate operator implementation. After a cylindrical lens integrates over the x-component, the final resulting state is measured in the y-coordinate. We test this using our basis and the Hadamard transform as example, and demonstrate the Deutsch–Jozsa algorithm with our encoding scheme, showing it can be used to query balanced and constant functions with average fidelity above 90%. In our scheme, we exploit the inherent parallelism offered by optical matrix-vector multiplication, allowing to prepare multiple identical spatial modes as a superposition and encode the oracle matrix as local phase transformations on each state.

We see that the matrix-vector multiplication can be decomposed into the elementwise product M ⊙ U and a subsequent summation of the columns. These operations can be achieved optically using spatial light modulators (SLMs) and a cylindrical lens (CL), as shown in Fig. 1(a), as demonstrated in Ref. 27. On one SLM, a uniform field can be modulated with a hologram that encodes (point-wise) the matrix U and is imaged onto a second SLM, which is encoded with the matrix M. The imaging system between the two SLMs enables the elementwise multiplication of the encoded matrices. The final operation is performed using a cylindrical lens that focuses the field in the direction that coincides with the columns of M, performing a column-wise summation—therefore completing the product to produce the elements v_j.

Remarkably, hidden in this approach is the fact that the transverse Cartesian coordinates are used as independent degrees of freedom that can be used to enact operations, analogous to non-interacting multi-particle systems in quantum computing. Next, we show how this idea can be used to emulate quantum computing with transverse fields.

Next, we show how this scheme can be used to perform a well-known quantum algorithm.

The Deutsch–Jozsa algorithm⁵⁰ provides an efficient quantum solution for determining whether a Boolean function, f(⋅), that takes binary values as inputs, is balanced (half of the inputs map to 1 and the other half to 0) or constant (all inputs either map to 0 or 1). To illustrate this, consider a collection of 2-bit input strings chosen from the set $B={00,01,10,11}$⁠; each element $k∈B$ can be used as an input to the function, i.e., F(k). A balanced function would return 0 for half the entries, e.g., k ∈ {00, 11}, and 1 for the other half, k ∈ {01, 10}, as shown in the first left panel of Fig. 2(a). The second example [top right panel of Fig. 2(a)] shows an alternative combination that can also represent a balanced function. For the constant case, all binary values map to a single output, i.e., {00, 01, 10, 11} → 0 or {00, 01, 10, 11} → 1 as shown in the bottom panel of Fig. 2(b) for two instances of such a function.

To query the nature of such a function, the quantum computer prepares a superposition of all possible binary values in $B$⁠, queries the function values simultaneously, and through interference it returns a result indicating whether the function is balanced or not. We adapt the algorithm to our encoding scheme using the lattice basis and our operator-state (matrix-vector) multiplication scheme.

In Fig. 3, we illustrate the setup used for our demonstration. To ensure a uniform field intensity across the transverse plane, we overfilled the first spatial light modulator (SLM₁) by magnifying a laser beam from a Helium Neon (HeNe) source with a wavelength of 633 nm using the telescope. On SLM₁, the input vector for the lattice of Gaussian modes was encoded as an array of displaced Gaussian beams, producing the superposition state, |Ψ⟩ ≡ |u₀⟩ from Eq. (2). To alter the coefficients of the lattice, each row of the matrix (U) mapping the vector was encoded using a digital hologram that encodes amplitudes and phases⁵¹ onto each of the Gaussian modes. Because the SLM encodes the vector state as a matrix, U, given any state |u⟩, the elements were copied into its rows following the procedure outlined in the theory.

The matrix generated on SLM₁ was then transferred to SLM₂ using an imaging system consisting of lenses L₁ and L₂. An aperture (AP) was placed between the two lenses to filter the modulated first order diffraction pattern from SLM₁. On SLM₂, we applied an additional phase encoding to the lattice/matrix elements, thereby completing the optical representation of the matrix-vector multiplication process. Following the encoding on SLM₂, the beam was passed through a cylindrical lens (CL) to undergo a one-dimensional Fourier transform. This operation focused the outcome of the matrix-vector multiplication into the central fringe (shown as an inset in Fig. 3), specifically at the zeroth order Fourier component in coordinate space. To extract this outcome accurately, we employed digital filtering, enabling precise capture of the resultant vector as shown in the inset of Fig. 3.

To characterize our system, we exploit the orthogonality of the Hadamard basis. In order to do this, we prepared the Hadamard basis states |u_j⟩ on SLM₁ and encoded the Hadamard transform on SLM₂. Thereafter, we propagated the field through the cylindrical lens and collected the resulting interference pattern on a CCD camera. The multiplication between the Hadamard basis states and the Hadamard gate enables us to clearly distinguish between different columns of the output vector state. Accordingly, we obtain orthogonal states (|v_j⟩) with distinct positions marked by bright lobes. We show a measured example of this for N = 2 dimensions in Fig. 4(a), where the first basis state, |u₀⟩, is encoded on the lattice and measured using the Hadamard gate. The outcome is a measured interference pattern, where the resulting filtered zero order produces a bright lobe in the first entry—for this specific example, the result corresponds to the state |v₀⟩. Summing and normalizing the intensity for each region of the output vector corresponded to extracting the output vector elements.

We repeated the process for all the basis states for a given dimension. In Fig. 4(b), we show the measured lobes for all the states in N = 2, confirming that the Hadamard basis states (|u₀⟩ and |u₁⟩) are mapped on to the basis states |v₀⟩ and |v₁⟩. From the measured lobes, we extracted the intensity from each component and produced a crosstalk matrix (CM). For N = 2, the crosstalk matrix resembles an identity matrix, showing that the measured basis modes are distinguishable and, therefore, confirms orthogonality between the measured vectors. To quantify this, we measured a fidelity of F = 0.99 ± 0.01, showing that the Hadamard matrix can map the basis elements |u_j⟩ onto the basis states |v_j⟩. This was calculated by normalizing each row and computing the average over the diagonal components.

Figures 4(c)–4(e) show the extension of this process to higher-dimensional matrices. For dimensions N = 4, N = 8, and N = 16, we obtain average fidelities of 0.96 ± 0.01, 0.94 ± 0.05, and 0.81 ± 0.02 multiplication, respectively. However, as the dimensions increase, we note that the fidelity decreases. We attribute this to the diffraction of Gaussian beams. Increasing the dimension or number of Gaussian modes for the same physical size of the matrix means decreasing the beam waist w₀ of each mode. As the beam waists get smaller, they expand at a faster rate during propagation, according to $w(z)=w01+λzπw022$⁠. Given that the cylindrical lens focuses light in one direction while neglecting the other, we expect the unfocused direction to continue expanding more rapidly as w₀ gets smaller, resulting in overlapping vector elements that are hard to distinguish. This diffraction effect impacts the fidelity of higher-dimensional computations.

Using our scheme, we demonstrate the Deutsch–Jozsa algorithm for a function that encodes information using N = 2, N = 4, and N = 8 basis elements, shown in Figs. 5(a)–5(c), respectively. We initialized our system in the superposition state |u₀⟩ by digitally displaying a hologram encoding the uniform lattice of Gaussian modes on SLM₁. The action of the oracle phase transformation and the subsequent Hadamard gate are encoded as a compound unitary transformation $HÛF$ on the second (SLM₂). The light beam is then sent through a cylindrical lens, which acts as a summing operator. The resulting interference pattern that was measured with the CCD camera and filtered zero order vector components are shown as insets in each figure.

From these images, we then plot the cross section of the normalized intensity against the vector element number. In this case, we can distinguish whether the function is balanced or constant by interrogating the presence or absence of light in the spatial position of the vector element. Using our basis vectors from the characterization, this means a constant function produces the vector |v₀⟩ and if the function is balanced |v_j≠0⟩. The instance of the oracle operator $(ÛF)$ used for each type of function is shown in each figure; this is shown as an inset in the intensity plots representing the measured components.

For each case, we accurately determined whether the function encoded in the unitary operation is constant or balanced. The average fidelity of these measurements exceeded 90% across all cases. In particular, for N = 2, 4, and 8, we achieved fidelities of 0.99 ± 0.01, 0.97 ± 0.01, and 0.93 ± 0.05, respectively.

In this study, we utilized optical matrix-vector multiplication to emulate quantum computing by treating the x–y components of a coherent field as a tensor product space, $HN⊗HN$⁠, derived from individual coordinates. The encoding basis was defined within this Hilbert space using a lattice of Gaussian modes, where the positions of the centers in the x–y coordinates established localized states. Vector states were encoded along the x-coordinate, while both coordinates were employed to encode the matrix operator. The output vector was measured along the y-coordinate after filtering out the zero-order component of the output interference pattern. We demonstrated the use of Hadamard basis elements and the Hadamard gate, achieving encoding and decoding fidelities as high as 95% for N = 8 dimensions and above 80% for N = 16.

Furthermore, our implementation of the Deutsch–Jozsa algorithm using optical vector matrix multipliers leverages the representation of quantum states as vectors and quantum operators as matrices that can be imprinted on our lattice basis. Optically, we prepare multiple identical Gaussian modes in parallel, with each mode representing an element of the states that are input to the function F(k). The high-dimensional nature of our encoding scheme enables the oracle function in our demonstration to be encoded as a diagonal operator.

However, we note spatial resources required to perform computations using structured light in higher dimensions are fundamentally limited by the diffraction limit. The minimum separation between distinguishable spatial modes is constrained by the wavelength of light and the numerical aperture of the optical system, leading to crosstalk between modes when attempting to scale up to higher dimensions. In our current setup, we have explored up to 16 distinguishable spatial modes, equivalent to a 4-qubit system. Spall’s vector-matrix demonstration reached up to ≈50 dimensions of encoding—in principle, every pixel of the SLM can be utilized as an encoding channel but would require much effort. Moreover, as the number of spatial modes increases, the complexity of addressing and manipulating these modes grows significantly. In quantum systems, the number of operations required for a computation scales polynomially with the number of qubits, whereas in classical systems, the need to individually resolve and address each spatial mode results in a linear scaling with the number of available modes. This ultimately limits the scalability of classical systems in performing quantum-like operations.

In addition, the preparation of Gaussian lattices and their unitary transformation matrices on an SLM typically involves NxN operations to generate the necessary hologram for the spot array, which can be computationally expensive and inefficient when scaling to higher dimensions. An alternative approach is the use of fan-out operations to generate the lattice of Gaussian modes.^52,53 Techniques such as cylindrical lenses⁵⁴ or diffractive optical elements⁵⁵ can directly produce the desired lattice by splitting or reshaping the input beam, reducing the need for pixel-by-pixel modulation on the SLM and significantly improving efficiency. Recent advances in trained diffractive optical elements provide another potential solution by allowing computationally designed optical masks to perform the necessary optical transformations.³⁵ This approach leverages inverse design techniques to optimize optical circuits, embedding high-dimensional transformations in complex media without requiring precise control over individual elements. Such masks streamline the process by reducing the need for real-time computation and re-calibration, ensuring scalability, and maintaining the accuracy of the transformations.

We acknowledge the proposal by Perez-Garcia et al.,⁴⁰ which employs classical light following Dragoman’s method. However, this approach necessitates multiple optical elements to prepare the superposition state. In contrast, our method efficiently prepares the state in a single step by projecting light onto a spatial light modulator (SLM) in a reprogrammable manner—the digital encoding allows for dynamic control and can, therefore, be adapted for any operation. Therefore, we anticipate that any transformation (operator/gate) can be encoded, including X-gate, Z-gate, etc., because the method allows any matrix operation to be encoded into the lattice.

The authors acknowledge Professor Andrew Forbes for his fruitful discussions and valuable advice. The authors also acknowledge the South African Quantum Initiative (SAQuti), the National Research Foundation (South Africa), and Optica for funding.

The authors have no conflicts to disclose.

Mwezi Koni: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Writing – original draft (equal). Hadrian Bezuidenhout: Investigation (supporting); Methodology (supporting); Writing – review & editing (supporting). Isaac Nape: Conceptualization (equal); Formal analysis (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.