A self-similar sine–cosine fractal architecture for multiport interferometers

1 Introduction

Photonic integrated circuits (PICs) have recently captured interest as a promising time- and energy-efficient platform for classical and quantum optical information processing. They have been used to accelerate tasks in signal processing [1–5], machine learning [6, 7], optimization [8], and quantum simulation [9–13]. Scaling these systems up in order to tackle real-world problems requires careful attention to issues such as the effect of analog component imperfections on performance, and the scaling of chip area with system size.

For instance, it has been shown that the test accuracy of optical neural networks (ONNs) based on Mach–Zehnder Interferometer (MZI) meshes [6] drops rapidly as soon as the constituent beamsplitters deviate from 50 to 50 splitting ratio by a couple of percent [14]. A variety of error-correction techniques have been proposed for the MZI-based platform—global optimization [15–20], local correction [21, 22], self-configuration [23–26], and hardware augmentation [27–29]. The behavior of many of these techniques can be better understood by considering an important insight derived in Ref. [21]—MZIs with imperfect beamsplitters implement only a subset of all the 2 × 2 unitary matrices that a perfect MZI can implement. This fact explains the observed imperfection-induced reduction in ONN performance since circuits composed of imperfect MZIs implement fewer functions than those with perfect MZIs [6, 14].

We show in this paper that the extent of reduction of the expressivity of a faulty MZI mesh depends strongly on its geometry and that a careful choice of mesh geometry can significantly soften the negative impact of hardware errors. We do so by introducing a novel self-similar MZI-mesh architecture based on the recursive sine–cosine unitary decomposition of Polcari [30] and demonstrating that it is more robust to MZI errors than the conventional Reck (triangular) [31] and Clements (rectangular) [32] mesh geometries. The recursive sine–cosine decomposition [30] is a generalization of the standard FFT decomposition of Fourier transform matrices [33] to arbitrary unitary matrices. We shall refer to MZI meshes constructed using this decomposition as sine–cosine fractal (SCF) meshes. Like the FFT mesh [34, 35], the SCF mesh has a recursive, self-similar structure; the FFT mesh can in fact be obtained from the SCF mesh by the mere pruning (omission) of certain columns of MZIs.

While SCF meshes have greater error robustness than other architectures,they can also be systematically shrunk in size for use in machine learning applications. In analogy with pruning in conventional neural networks [36, 37], we introduce a systematic mesh pruning scheme that interpolates between the simple FFT and the full sine–cosine fractal, and numerically demonstrate that ONNs composed of pruned meshes still achieve excellent performance at benchmark learning tasks.

The paper is organized as follows: Section 2 introduces and discusses the sine–cosine fractal architecture. Section 3 contains analytical and numerical results on the expressivity of SCF meshes in the presence of beamsplitter errors. Section 4 reports the performance of ONNs constructed from both complete and pruned SCF meshes, and Section 5 concludes the paper with a further discussion of the scope and impact of our work.

2 The sine–cosine fractal architecture

In order to construct a photonic circuit that implements a given unitary matrix U, one first decomposes U into a product of 2 × 2 unitary matrices and diagonal phase shifts. MZIs are used to implement the 2 × 2 unitary matrices in the hardware; the transfer function of an MZI with two phase-shifters θ, ϕ (Figure 1(a)) is given by:

The arrangement of MZIs in the circuit, that is, the geometry of the MZI mesh, determines the order of appearance of the corresponding 2 × 2 unitaries in the decomposition. Figure 1(a) depicts a mesh that implements an 8 × 8 matrix via the Clements decomposition [32]. The Reck [31] and balanced binary tree [38] are other important decompositions that have been used to construct unitary meshes.

Figure 1:

Illustration of optical mesh designs. (a) 8 × 8 Clements mesh. (b) First step of the block decomposition of the mesh. This results in 4 × 4 quadrants. (c) Further decomposition of the quadrants, which results in the SCF mesh.

This paper proposes a new mesh architecture based on the sine–cosine decomposition, a block diagonalization of unitary matrices, which on an N × N matrix U consists of partitioning U into four N/2 × N/2 blocks and performing the singular value decomposition (SVD) on each block. The unitarity of U imposes special constraints that force the blocks to share singular vectors. Polcari [30] shows that the block-wise SVD of U yields:

where U₁₁, U₁₂, U₂₁, U₂₂ are unitary matrices and D₁₁, D₁₂, D₂₁, D₂₂ are diagonal matrices encoding the singular vectors and values, respectively. Unitarity constrains the D _ij to take the following form:

with (Θ, Φ) representing diagonal matrices that encode phase shifts. Figure 1(b) depicts Eq. (2) graphically for the 8 × 8 case—the 4 × 4 unitary matrices U _ij are implemented by Clements meshes while the diagonal matrices D _ij are implemented by MZIs in the center that couple the four unitary blocks. One can actually go further and perform the block-wise SVD again on each of the U _ij sub-blocks to obtain the mesh of Figure 1(c); the 2 × 2 unitaries obtained from the 4 × 4 unitaries are now directly implemented by MZIs. Because of its self-similar structure, we denote this geometry the sine–cosine fractal (SCF) mesh. In the general case, an SCF mesh can be constructed from any radix-2 (N = 2 ⁿ ) matrix: one recursively performs blockwise SVDs on each unitary matrix of size greater than 2 until the full decomposition consists only of 2 × 2 matrices that connect different modes.

Like the Clements and other conventional architectures, the SCF mesh has a depth that scales as O(N), has N² degrees of freedom, and is universal, i.e. it can represent the entire unitary group. However, the SCF mesh also possesses tunable crossings that couple non-neighboring waveguides (Figure 1(c)), i.e. crossings of stride s > 1, interleaved between conventional crossings with s = 1. This is in contrast to conventional architectures, where all MZIs have unit stride.

3 Error correction and matrix fidelity

The long-stride crossings give the SCF mesh greater robustness to hardware imperfections. To see how, consider the problem of realizing a target at high fidelity on imperfect hardware. As mentioned previously, MZI meshes with perfect beamsplitters can implement any unitary matrix; the introduction of faults, however, reduces the expressivity of the mesh and consequently the fraction of matrices that are implementable drops below unity. In this section, we show that sine–cosine fractal meshes can perfectly implement a greater fraction of random matrices than Clements meshes can for the same beamsplitter error level.

3.1 Distribution of mesh phase-shifts for Haar-random matrices

A specific setting of phase shifts (θ, ϕ) is required to make a mesh implement a given target unitary matrix. Drawing the target matrix from the Haar (uniform) distribution [39] induces a distribution P(θ) over the phase-shifts. Russell et al. [40] show that the phase-shift distribution of the nth MZI (according to any indexing scheme) in either the Reck or Clements meshes is given by P n ( θ ) = k n ⁡ sin ( θ / 2 ) cos ( θ / 2 ) 2 k n − 1 , where k _n ∈ {1, …, N − 1}, called the rank of the nth MZI, is a function of the physical position of the MZI within the mesh. There are (N − k) MZIs of rank k [40]. For larger meshes, the rank of an average MZI increases, and this results in P(θ), which is an average of the P _n (θ) over all the MZIs, clustering around the cross state θ = 0 (top row of Figure 2). The phase shift distribution for both the Clements and Reck meshes are identical, with P(θ) clustering near θ = 0 [40].

Figure 2:

Distribution of crossing angles for the Clements, Reck and SCF meshes as a function of mesh size and stride, respectively. The crossing angles of the Clements mesh and the Reck mesh increasingly clusters near θ = 0 for larger meshes. Crossing angles for the SCF mesh become increasingly uniformly distributed as the stride increases.

However, unlike the Reck and Clements meshes, the SCF mesh is configured from a top–down block decomposition of the matrix. For a given matrix U, sampled over the Haar measure, the singular-vector matrices U _ij , (i, j) ∈ {1, 2} in Eq. (2) are also Haar-random and independent of each other. As a result, the distribution P _n (θ) depends only on the stride s _n of the nth MZI, not on its location in the mesh. The bottom row of Figure 2 shows the distribution of angles for MZIs of stride s = 1, 4, 32 for mesh size N = 64. The majority of MZIs have unit stride and P _n (θ) ∝ sin(θ). As the stride increases, this distribution begins to resemble a uniform distribution. We find that P _n (θ) for stride s _n can be fit well by the normalized finite Fourier series of a constant given by:

(4) P n ( θ ) ≈ 1 2 ∑ q = 1 s n ( 2 q − 1 ) − 2 ∑ q = 1 s n sin ( 2 q − 1 ) θ 2 q − 1

3.2 Error in implementing Haar-random matrices in the presence of MZI errors

The deviations from 50:50 of the two constituent beamsplitters of fabricated MZIs are captured by the phase angles (α, β). These splitter errors perturb the MZI transfer matrix as follows:

(5) T ′ ( θ ′ , ϕ ′ , α , β ) = cos π 4 + β i sin π 4 + β i sin π 4 + β cos π 4 + β e i θ ′ 0 0 1 × cos π 4 + α i sin π 4 + α i sin π 4 + α cos π 4 + α e i ϕ ′ 0 0 1

Bandyopadhyay et al. [21] show that it is always possible to choose phase-shifts θ′, ϕ′ for a faulty MZI with errors (α, β) such that it implements the transfer matrix of an ideal MZI T(θ, ϕ) as long as:

(6) 2 | α + β | ︸ θ min ≤ θ ≤ π − 2 | α − β | ︸ θ max

If θ is outside this range, the faulty MZI cannot exactly emulate the ideal MZI and the faulty mesh transfer function deviates from that of the ideal mesh which implements the target matrix. Alternate approaches to correct for these errors have been proposed in refs. [15, 27, 29, 41].

For a given target unitary U, we quantify the deviation of the faulty mesh using the Frobenius norm  E = ‖ Δ U ‖ / N which computes the average relative error per matrix element. This quantity is then averaged over both the choice of target unitary (from the Haar distribution), and the distribution of MZI splitter errors, which are assumed to be independent Gaussians ( α , β ∼ N ( 0 , σ ) ) . In the case of correlated errors, the effect of most correlations vanishes over the Haar measure as proven in Ref. [24].

In the absence of error-correction techniques, all the MZIs of the ideal mesh are implemented incorrectly by the corresponding MZIs in the faulty mesh and the error Frobenius norm  E 0 (for both Clements and SCF meshes) has been shown [24] to scale linearly with the MZI error σ as E 0 = 2 N σ . When one uses the error-correction techniques of Refs. [21, 24, 25], only the MZIs in the ideal mesh that do not satisfy Eq. (6) are implemented incorrectly by the corresponding MZIs in the faulty mesh—only those MZIs contribute to the error Frobenius norm. Using n to denote the MZI location in the mesh as before, the average “corrected” error, E c (assuming uncorrelated errors, see [27, Supp. Sec. 1]), is computed by integrating over the probability that each MZI does not satisfy Eq. (6):

(7) ( E c ) 2 = 1 2 N ∑ n ∫ 0 θ min P n ( θ ) ( θ − θ min ) 2 d θ   + ∫ θ max π P n ( θ ) ( θ − θ max ) 2 d θ

Eq. (7) indicates that the effectiveness of error-correction, measured by E c , is strongly dependent on the distribution P _n (θ). For the Clements mesh, Ref. [24] proves that E c ( c l e m ) = 2 3 N σ 2 , which is a quadratic improvement over E 0 .

Since the integral is over angles close to either θ = 0 or θ = π, the two terms in Eq. (7) are estimated for the sine–cosine fractal mesh by Taylor-expanding Eq. (4) to first order about θ = 0, π, respectively. The result is:

(8) E c ( s c f ) = 4 π N ⁡ log 2 ( N ) σ 2 .

The ratio of corrected errors for both meshes is:

(9) E c ( c l e m ) E c ( s c f ) = π 2 24 N log 2 ( N )

which is greater than 1 for all but very small N.

We performed numerical experiments on meshes up to size N = 1024 to validate the above expressions (Figure 3). The observed uncorrected error E 0 of both meshes (Figure 3(a)) scales as σ while error-correction on both meshes improves this scaling to σ². For meshes of size 256 × 256, error correction shows over an order of magnitude improvement in matrix error Frobenius norm over the uncorrected case when σ ≤ 0.01, with the Sine–Cosine Fractal mesh performing better than the Clements mesh. Figure 3(b) illustrates the distribution of post-correction error Frobenius norm as a function of mesh size N for a fixed σ = 0.02 (which is a typical value for directional couplers under wafer-scale process variation). While the factor N log 2 ( N ) is modest for small meshes (N < 64), it clearly opens up a significant accuracy gap in Figure 3(b) in the large scale regime (N > 1024).

Figure 3:

Comparision of matrix error scaling for the SCF and Clements meshes. (a) Scaling of matrix error with process variation σ. (b) Scaling of matrix error with mesh size, showing the advantage of the SCF mesh over the Clements mesh for larger mesh sizes.

3.3 Fraction of Haar-random matrices those are exactly realizable in the presence of MZI errors

The fraction of the unitary group U(N) that can be realized by an imperfect mesh is equal to the probability that, under the Haar measure [40], all target splitting angles θ are realizable. For convenience, we derive a quantity called the coverage, cov(N) [24], from this probability:

(10) cov ( N ) = ∏ n 1 − P n ( θ < θ min ) − P n ( θ > θ max ) ≈ exp − ∑ n P n ( θ < θ min ) + P n ( θ > θ max )

For a Clements mesh, Ref. [24] proves that cov ( clem ) ( N ) = exp − N 3 σ 2 3 . Like the matrix error in the previous subsection, the coverage of an SCF mesh is computed using Taylor series expansions. The result is:

(11) cov ( scf ) ( N ) = exp − 8 N 2 ⁡ log 2 ( N ) π 2 σ 2

which is greater than cov^(clem)(N) for the same σ for all but very small N.

This coverage, however, expresses the ability of an imperfect mesh to realize random matrices under the Haar measure. It is often the case that such meshes will be used to realize matrices that are very commonly used in quantum computation, signal processing or machine learning. In Figure 4(a), we compare the ability of the SCF and Clements meshes to realize Hadamard matrices and discrete-Fourier transform (DFT) matrices of size N = 128. These examples confirm that many common matrices possess structure that does not follow the statistics of Haar random unitaries. For the Hadamard matrix, the Clements decomposition yields a mesh where 50% of all MZIs are programmed to the crossing state θ = 0 (Figure 4(b)); as a result, the matrix fidelity is much worse than for a typical Haar-random instance, and scales as E c ∝ σ as opposed to σ². The SCF mesh, by contrast, has a more uniform distribution of splitter angles, recovering E c ∝ σ 2 scaling with overall matrix errors about an order of magnitude smaller.

Figure 4:

Realization of Hadamard and DFT matrices with Clements and SCF meshes. (a) Post-correction error E c for matrix size N = 128. (b) MZI splitter angles θ for meshes of size N = 64 realizing Hadamard and DFT matrices.

For the DFT matrix, the situation is reversed. Here, the Clements mesh shows a more uniform distribution of angles, while the SCF mesh has a large number of MZIs with extreme splitting ratios θ ∈ {0, π}. Therefore, the SCF mesh achieves significantly worse matrix fidelity when straightforwardly implementing the DFT. The reason is connected with the highly non-uniform distribution of block-singular values (i.e. D _ij in Eq. (3)) in the standard DFT. We can circumvent this problem using the technique of port allocation [22]. A random permutation of the input and output ports yields a more even distribution of singular values, which reduces the number of MZIs with θ ∈ {0, π}. The resulting matrix fidelity is comparable to that of the Clements mesh (green curve in Figure 4(a), right).

4 Use in optical neural networks

In this section, we study the performance of an optical neural network (ONN) built from sine–cosine fractal meshes. We also propose and evaluate a pruning scheme for these meshes that allows areal footprint reduction while maintaining test performance. The neural network configuration is similar to those studied in Refs. [6, 14, 26]—each neural net layer is implemented by an SCF mesh connected to an electro-optic nonlinearity [42] (Figure 6(a)). All our networks had two layers. Our simulations used the meshes [43] package, and results are presented for the MNIST [44] image classification task. The preprocessing of the images involved low-pass filtering and was identical to the procedure adopted in Refs. [21, 26]. The standard cross entropy loss and the Adam optimizer were used for training.

ONNs trained with SCF meshes achieved accuracies that matched the Clements mesh [21, 24] of ∼95%–96% for small meshes (N = 64) and ∼97% for larger meshes (N = 256). Next, we simulated the effect of MZI errors on the trained SCF mesh neural net; Figure 6(b) shows the median classification accuracy of 10 independently trained networks as a function of splitter errors. SCF networks of size N = 64 yielded ∼95% test accuracy while those of size N = 256 reached ∼97%. The presence of MZI errors rapidly degrades the performance of the network, with accuracy dropping to below 90% with splitter variation as low as ∼2%. The use of hardware error correction, however, extends this cutoff to greater than 12% even for bigger meshes, which is well above present-day process error [45] and larger than the corresponding cutoff for Clements meshes (which is 6% [21]).

4.1 Weight pruning

The number of columns of MZIs with stride s is N 2 s in the Sine–Cosine Fractal mesh while the standard FFT mesh contains a single column of each stride. We introduce a pruning scheme that interpolates between these extremes by introducing a fractal dimension D ∈ [ 1,2 ] —in a partially pruned mesh, the number of columns of stride s is N 2 s D − 1 . Setting D = 1 and D = 2 yields the FFT and SCF meshes, respectively. Controlling D allows us to tune the number of degrees of freedom and reduce the areal footprint of the device while ensuring full connectivity. Figure 5 illustrates partially pruned 8 × 8 SCF meshes for different values of D . The depth of the pruned network (approximated to leading order in N) scales as N D − 1 2 D − 1 − 1 while the number of MZIs is approximately N D 2 2 D − 1 . To construct a pruned mesh, the column with the largest stride is chosen initially and is interleaved with columns of the next smaller stride. Each of these columns is then interleaved with the following smaller stride from left to right until s = 1 is reached. When the maximum number of columns is reached for a given stride, the interleaving is terminated and the process is continued with a column of the next smaller stride. This ensures that no two consecutive columns have the same stride since any such columns would collapse into a single column.

Figure 5:

The SCF mesh at different stages of pruning constructed by interleaving strides from the largest successively to the smallest, from left to right with the number of columns limited by N 2 s D − 1 : (a) A maximally pruned mesh, which is the standard FFT mesh, D = 1.0 . (b) D = 1.5 (c) D = 1.75 .

Figure 6(c) illustrates the results of training networks with pruned meshes of ideal MZIs for different D . Increasing D (decreasing the amount of pruning, or increasing the size of the network) increases the classification accuracy as one would expect. Interestingly, a maximally pruned (that is, a standard FFT) 64 × 64 2-layer ONN already achieves 95% accuracy, which is commensurate with the performance of present-day DNN accelerators [46].

Figure 6:

Performance of optical neural networks (ONNs) constructed using SCF meshes. (a) 2-Layer deep neural network architecture for MNIST classification where the blocks U^(0,1) represent unitary transforms by the sine–cosine fractal mesh. (b) Simulated classification accuracy as a function of mesh size (N = 64 and 256) and MZI error σ. (c) Classification accuracy as a function of degree of pruning D ∈ [ 1,2 ] and mesh size (N = 16 and 64 and 256). (d) and (e) Classification accuracy as a function of the degree of pruning D and MZI error, trained on maximally faulty mesh of size N = 64 and 256.

5 Conclusions

This work presented a novel architecture for multiport interferometers based on the sine–cosine decomposition of unitary matrices. The proposed scheme is self-similar, and therefore modular. As a result, this architecture is ideal to construct multi-chiplet modules for large-scale devices that are typically limited by device yield. We showed that SCF meshes show improved scaling under error correction when compared to conventional multiport interferometers. Finally, this design allows for systematic re-wiring of MZI layers, which is an efficient way of reducing the areal footprint of the mesh while maintaining full connectivity.

The proposed design has multiple advantages over the traditional architectures for multiport interferometers. Due to the uniform distribution of coupling angles (which is in contrast to conventional mesh architectures where crossing angles are clustered around the cross state), error correction techniques are far more effective in the case of SCF meshes. Truncated meshes used in ONNs can be trained to perform on par with the performances of present-day DNN accelerators. While the benefits of modularity and stronger robustness to error are small for present-day mesh size, the scaling gap of O(N) versus O ( N ⁡ log 2 ( N ) ) significantly impacts the performance for larger meshes. This reduced scaling with N implies that SCF meshes are more sensitive to improvements in foundry process. Reducing σ by 2× corresponds to a 4× increase in the maximum mesh size for Reck and Clements owing to the O(Nσ²) scaling; for the butterfly fractal, the corresponding increase would be 16×.

A potential drawback of this scheme is the presence of multiple large stride crossings with non-zero crosstalk. The loss/crosstalk introduced by these crossings can be minimized by the following methods that are enabled by the modularity of the SCF architecture:

1. Incomplete Decomposition: To minimize the number of crossings, the block decomposition of the mesh can be terminated such that the smallest block size N_blk > 2. Each block will then take the form of a standard Clements geometry with no intra-block crossings.

2. Out-of-plane Crossings: Inter-block crossings (s > 1 among arbitrary block sizes) can be implemented using waveguide escalators and out-of-plane crossings. Crossings of different stride can be fabricated into multiple layers of the photonic circuit and can be connected with inter-planar couplers. The structure of the SCF mesh ensures that each spatial mode would travel through the same number of inter-planar couplers, which results in losses that are uniform. This would only impact the power scaling of the output vector, and not the matrix realized by the mesh. Recent work on the construction of these couplers has shown losses as low as 0.05 ± 0.02 dB or approximately 98.7% per coupler [48]. Scaling up the input power by a factor of 1/0.987^(N−2) (up to heating limitations posed by the hardware) for a mesh of size N would allow us to account for losses due to the couplers. Out-of-plane crossings have also been shown to have much lower crosstalk than in-plane crossings [49, 50]. Crosstalk on the order of −50 dB for input powers of 1 mW has been achieved in non-planar topologies, and can be fabricated scalably for such an integrated photonic mesh [51, 52]. Since the crosstalk due to induced fields would decrease exponentially with increasing interplanar distance, it would be easily possible to further reduce the crosstalk. In the case of a multi-chip module, each chiplet will be connected by a set of waveguide crossings with a large stride. These crossings can be fabricated using present-day lithography and laser-writing techniques in either polymer [53, 54] or glass [55, 56], and could utilize “hockey-stick” escalator couplers to reduce the alignment tolerances for each chiplet [57].

For practical use as energy efficient deep learning accelerators, nanophotonic circuits will need to be scaled to reach the sub-fJ/MAC energy target. Present-day foundry processes face immense difficulties in scaling conventional multiport interferometers to the large size required to achieve this energy efficiency target. Furthermore, these meshes suffer unacceptably high errors E c ( c l e m ) ≳ 0.2 even in the presence of error correction [21, 24, 25]. This brings us to the regime below the recommended 4–8 bits of precision that is usually targeted for DNN accelerators [58–60]. The modularity of the SCF mesh, and its improved tolerance to error, open up the regime of large-scale programmable photonics, making the SCF mesh frontrunning candidate architecture for future systems.