\section{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 \cite{annoni2017unscrambling, ribeiro2016demonstration, milanizadeh2019manipulating, zhuang2015programmable, notaros2017programmable}, machine learning \cite{NN_arch_1, opt_ml}, optimization \cite{prabhu2020ising}, and quantum simulation \cite{harris2017quantum, wang2018multidimensional, qiang2018large, sparrow2018simulating, carolan2015universal}. 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 \cite{NN_arch_1} drops rapidly as soon as the constituent beamsplitters deviate from 50-50 splitting ratio by a couple of percent \cite{NN_arch_2}. A variety of error-correction techniques have been proposed for the MZI-based platform\textemdash global optimization \cite{imperfect_mesh_1, imperfect_mesh_2, lopez2019programmable, lopez2020auto, perez2020multipurpose, imperfect_mesh_3}, local correction \cite{hec_1, imperfect_mesh_4}, self-configuration \cite{imperfect_mesh_5, self_con_1, self_con_2, NN_arch_3}, and hardware augmentation \cite{Inf, suzuki2015ultra, miller2015perfect}. The behavior of many of these techniques can be better understood by considering an important insight derived in Ref.~\cite{hec_1}\textemdash MZIs with imperfect beamsplitters implement only a subset of all the $2\times2$ 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 \cite{NN_arch_1, NN_arch_2}.

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 \cite{butterfly_decomp} and demonstrating that it is more robust to MZI errors than the conventional \texttt{RECK} (triangular) \cite{reck} and \texttt{CLEMENTS} (rectangular) \cite{clements} mesh geometries. The recursive Sine-Cosine decomposition \cite{butterfly_decomp} is a generalization of the standard FFT decomposition of Fourier Transform matrices~\cite{butterfly_og} to arbitrary unitary matrices. We shall refer to MZI meshes constructed using this decomposition as Sine-Cosine Fractal (SCF) meshes. Like the FFT mesh~\cite{FFT_analog_1, FFT_analog_2}, 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~\cite{conventional_prune_1, conventional_prune_2}, 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~\ref{sec:2} introduces and discusses the Sine-Cosine Fractal architecture. Section~\ref{sec:3} contains analytical and numerical results on the expressivity of SCF meshes in the presence of beamsplitter errors. Section~\ref{sec:4} reports the performance of ONNs constructed from both complete and pruned SCF meshes, and Section~\ref{sec:5} concludes the paper with a further discussion of the scope and impact of our work.

\section{\label{sec: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\times 2$ unitary matrices and diagonal phase shifts. MZIs are used to implement the $2\times 2$ unitary matrices in the hardware; the transfer function of an MZI with two phase-shifters $\theta,\phi$ (Fig.~\ref{fig_1}(a)) is given by:
\begin{align}
	T(\theta,\phi) & = \frac{1}{2} \begin{bmatrix} 1 & i \\ i & 1 \end{bmatrix} \begin{bmatrix} e^{i\theta} & 0 \\ 0 & 1 \end{bmatrix}
				  	    \begin{bmatrix} 1 & i \\ i & 1 \end{bmatrix} \begin{bmatrix} e^{i\phi} & 0 \\ 0 & 1 \end{bmatrix}
					    \nonumber \\
	  & = i e^{i\theta/2} \begin{bmatrix} e^{i\phi} \sin(\theta/2) & \cos(\theta/2) \\ e^{i\phi} \cos(\theta/2) & \sin(\theta/2) \end{bmatrix} \label{eq:t}
\end{align}
The arrangement of MZIs in the circuit, that is, the geometry of the MZI mesh, determines the order of appearance of the corresponding $2\times 2$ unitaries in the decomposition. Fig.~\ref{fig_1}(a) depicts a mesh that implements an $8\times 8$ matrix via the \texttt{CLEMENTS} decomposition~\cite{clements}. The \texttt{RECK}~\cite{reck} and balanced binary tree~\cite{binary_tree_decomp} are other important decompositions that have been used to construct unitary meshes. 

This paper proposes a new mesh architecture based on the Sine-Cosine decomposition, a block diagonalization of unitary matrices, which on an $N\times N$ matrix $U$ consists of partitioning $U$ into four $N/2\times 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 \cite{butterfly_decomp} shows that the block-wise SVD of $U$ yields:
\begin{equation}
    U = \begin{bmatrix} 
    U_{12} &0 \\ 
      0    & U_{22}
    \end{bmatrix}
    \begin{bmatrix}
    D_{11} & D_{12} \\ 
    D_{21} & D_{22}
    \end{bmatrix}
    \begin{bmatrix}
    U_{11} &0 \\
    0 & U_{21}
    \end{bmatrix}
    \label{eq:udiv}
\end{equation}
where $U_{11}, U_{12}, U_{21}, U_{22}$ are unitary matrices and $D_{11}, D_{12}, D_{21}, D_{22}$ are diagonal matrices encoding the singular vectors and values, respectively.  Unitarity constrains the $D_{ij}$ to take the following form: 
\begin{equation}
	\begin{bmatrix} D_{11} & D_{12} \\ D_{21} & D_{22} \end{bmatrix} = i e^{i\Theta/2} \begin{bmatrix} e^{i\Phi}\sin(\Theta/2) & \cos(\Theta/2) \\ e^{i\Phi}\cos(\Theta/2) & -\sin(\Theta/2) \end{bmatrix} \label{eq:dij}
\end{equation}
with $(\Theta, \Phi)$ representing diagonal matrices that encode phase shifts. Fig.~\ref{fig_1}(b) depicts Eq.~\eqref{eq:udiv} graphically for the $8\times8$ case\textemdash the $4\times4$ 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}$ subblocks to obtain the mesh of Fig.~\ref{fig_1}(c); the $2\times2$ unitaries obtained from the $4\times4$ 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^n$) matrix: one recursively performs blockwise SVDs on each unitary matrix of size greater than $2$ until the full decomposition consists only of $2\times2$ matrices that connect different modes.      

\begin{figure}
    \centering
    \includegraphics[width = \columnwidth]{Figure_pdfs/Fig_1_final.pdf}
    \caption{\textbf{(a)} $8 \times 8$ \texttt{CLEMENTS} mesh. \textbf{(b)} First step of the block decomposition of the mesh. This results in $4 \times 4$ quadrants. \textbf{(c)} Further decomposition of the quadrants, which results in the SCF mesh. }
    \label{fig_1}
\end{figure}

Like the \texttt{CLEMENTS} and other conventional architectures, the SCF mesh has a depth that scales as $O(N)$, has $N^{2}$ 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 (Fig.~\ref{fig_1}(c)), i.e.\ crossings of {\it stride} $s > 1$, interleaved between conventional crossings with $s = 1$.  This is in contrast to conventional architectures, where all MZIs have unit stride.

\section{\label{sec: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; however, 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.  

\subsection{Distribution of mesh phase-shifts for Haar-random matrices}

A specific setting of phase shifts $(\theta, \phi)$ is required to make a mesh implement a given target unitary matrix. Drawing the target matrix from the Haar (uniform) distribution \cite{Haar_sampling} induces a distribution $P(\theta)$ over the phase-shifts. Russell et al. \cite{MZI_rank} show that the phase-shift distribution of the $n$-th MZI (according to any indexing scheme) in either the \texttt{RECK} or \texttt{CLEMENTS} meshes is given by $P_n(\theta) = k_n \mathrm{sin}(\theta/2) \mathrm{cos}(\theta/2)^{2k_n - 1}$, where $k_n \in \{1, \hdots, N - 1 \}$, called the \textit{rank} of the $n$-th MZI, is a function of the physical position of the MZI within the mesh. There are $(N - k)$ MZIs of rank $k$~\cite{MZI_rank}. For larger meshes, the rank of an average MZI increases, and this results in $P(\theta)$, which is an average of the $P_n(\theta)$ over all the MZIs, clustering around the cross state $\theta = 0$ (top row of Fig.~\ref{fig_2}). 

However, unlike the \texttt{RECK} and \texttt{CLEMENTS} meshes, the SCF mesh is configured from a top-down block decompsition of the matrix.  For a given matrix $U$, sampled over the Haar measure, the singular-vector matrices $U_{ij}, (i,j) \in \{1, 2\}$ in Eq.~(\ref{eq:udiv}) are also Haar-random and independent of each other. As a result, the distribution $P_n(\theta)$ depends only on the stride $s_n$ of the $n$-th MZI, not on its location in the mesh. The bottom row of Fig.~\ref{fig_2} shows the distribution of angles for MZIs of different stride.  The majority of MZIs have unit stride and $P_n(\theta) \propto \sin(\theta)$. As the stride increases, this distribution begins to resemble a uniform distribution. We find that $P_n(\theta)$ for stride $s_n$ can be fit well by the normalized finite Fourier series of a constant given by:
\begin{equation}
    	P_n(\theta) \approx \frac{1}{2\sum_{q=1}^{s_n}{(2q-1)^{-2}}} \sum_{q=1}^{s_n}{\frac{\sin\bigl((2q-1)\theta\bigr)}{2q-1}}\label{eq:scffourier}
\end{equation}

\subsection{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 ($\alpha, \beta$). These splitter errors perturb the MZI transfer matrix as follows: 
\begin{align}
	T'(\theta',\phi',\alpha,\beta) & = 
	\begin{bmatrix} \cos(\tfrac\pi4\!+\!\beta) & i\sin(\tfrac\pi4\!+\!\beta) \\ i\sin(\tfrac\pi4\!+\!\beta) & \cos(\tfrac\pi4\!+\!\beta) \end{bmatrix} \begin{bmatrix} e^{i\theta'} & 0 \\ 0 & 1 \end{bmatrix} \nonumber \\
	& \ \ \ \ \times 	\begin{bmatrix} \cos(\tfrac\pi4\!+\!\alpha) & i\sin(\tfrac\pi4\!+\!\alpha) \\ i\sin(\tfrac\pi4\!+\!\alpha) & \cos(\tfrac\pi4\!+\!\alpha) \end{bmatrix} \begin{bmatrix} e^{i\phi'} & 0 \\ 0 & 1 \end{bmatrix} \label{eq:dt}
\end{align}
Bandyopadhyay et al. \cite{hec_1} show that it is always possible to choose phase-shifts $\theta',\phi'$ for a faulty MZI with errors ($\alpha, \beta$) such that it implements the transfer matrix of an ideal MZI $T(\theta,\phi)$ as long as:
\begin{equation}
    \underbrace{2|\alpha+\beta|}_{\theta_{\rm min}} \leq \theta \leq \underbrace{\pi - 2|\alpha-\beta|}_{\theta_{\rm max}} \label{eq:thetarange}
\end{equation}
If $\theta$ 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.

\begin{figure}
    \centering
    \includegraphics[width = \columnwidth]{Figure_pdfs/fig_2_final.pdf}
    \caption{Distribution of crossing angles for the \texttt{CLEMENTS} and SCF meshes as a function of mesh size and stride respectively. The crossing angles of the \texttt{CLEMENTS} mesh increasingly clusters near $\theta = 0$ for larger meshes. Crossing angles for the SCF mesh become increasingly uniformly distributed as the stride increases.}
    \label{fig_2}
\end{figure}

For a given target unitary $U$, we quantify the deviation of the faulty mesh using the Frobenius norm $\mathcal{E} = \lVert\Delta U\rVert/\sqrt{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 ($\alpha, \beta \sim \mathcal{N}(0, \sigma)$). In the case of correlated errors, the effect of most correlations vanish over the Haar measure as proven in Ref.~\cite{self_con_1}. 

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 $\mathcal{E}_{0}$ (for both \texttt{CLEMENTS} and SCF meshes) has been shown~\cite{self_con_1} to scale linearly with the MZI error $\sigma$ as $\mathcal{E}_{0} = \sqrt{2N}\sigma$. When one uses the error-correction techniques of Refs.~\cite{hec_1, self_con_1, self_con_2}, only the MZIs in the ideal mesh that do not satisfy Eq.~\eqref{eq:thetarange} are implemented incorrectly by the corresponding MZIs in the faulty mesh\textemdash 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, $\mathcal{E}_{c}$ (assuming uncorrelated errors, see \cite[Supp.~Sec.~1]{Inf}), is computed by integrating over the probability that each MZI does not satisfy Eq.~\eqref{eq:thetarange}:
\begin{align}
    	(\mathcal{E}_c)^2 & = \frac{1}{2N}\sum_n \Bigl[\int_{0}^{\theta_{\rm min}} P_n(\theta)(\theta-\theta_{\rm min})^2\mathrm{d}\theta \nonumber \\
    		& \qquad\qquad\ \  + \int_{\theta_{\rm max}}^{\pi} P_n(\theta)(\theta-\theta_{\rm max})^2\mathrm{d}\theta\Bigr] 
    \label{eq:ec}
\end{align}
Eq.~\eqref{eq:ec} indicates that the effectiveness of error-correction, measured by $\mathcal{E}_c$, is strongly dependent on the distribution $P_{n}(\theta)$. For the \texttt{CLEMENTS} mesh, Ref.~\cite{self_con_1} proves that $\mathcal{E}_{c}^{\mathrm{(clem)}} = \sqrt{\frac{2}{3}} N\sigma^{2}$, which is a quadratic improvement over $\mathcal{E}_0$. 

Since the integral is over angles close to either $\theta=0$ or $\theta=\pi$, the two terms in Eq.~\eqref{eq:ec} are estimated for the Sine-Cosine Fractal mesh by Taylor-expanding Eq.~\eqref{eq:scffourier} to first order about $\theta=0,\pi$ respectively. The result is:
\begin{equation}
    \mathcal{E}_c^{\mathrm{(scf)}} = \frac{4}{\pi} \sqrt{N\log_2(N)}\, \sigma^2. \label{eq:ecft}
\end{equation}
The ratio of corrected errors for both meshes is:
\begin{equation}
    \frac{\mathcal{E}_c^{\rm (clem)}}{\mathcal{E}_c^{\rm (scf)}} = \sqrt{\frac{\pi^2}{24} \frac{N}{\log_2(N)}}
\end{equation}
which is greater than $1$ for all but very small $N$.

\begin{figure}[t]
    \centering
    \includegraphics[width = \columnwidth]{Figure_pdfs/Fig_3_final.pdf}
    \caption{\textbf{(a)} Scaling of matrix error with process variation $\sigma$. \textbf{(b)} Scaling of matrix error with mesh size, showing the advantage of the SCF mesh over the \texttt{CLEMENTS} mesh for larger mesh sizes.}
    \label{fig_3}
\end{figure}

We performed numerical experiments on meshes up to size $N = 1024$ to validate the above expressions (Fig.~\ref{fig_3}). The observed uncorrected error $\mathcal{E}_0$ of both meshes (Fig.~\ref{fig_3}(a)) scales as $\sigma$ while error-correction on both meshes improves this scaling to $\sigma^{2}$. For meshes of size $256 \times 256$, error correction shows over an order of magnitude improvement in matrix error Frobenius norm over the uncorrected case when $\sigma \leq 0.01$, with the Sine-Cosine Fractal mesh performing better than the \texttt{CLEMENTS} mesh. Fig.~\ref{fig_3}(b) illustrates the distribution of post-correction error Frobenius norm as a function of mesh size $N$ for a fixed $\sigma = 0.02$ (which is a typical value for directional couplers under wafer-scale process variation). While the factor $\sqrt{\frac{N}{\mathrm{log}_{2}(N)}}$ is modest for small meshes ($N < 64$), it clearly opens up a significant accuracy gap in Fig.~\ref{fig_3}(b) in the large scale regime ($N > 1024$). 

\subsection{Fraction of Haar-random matrices that 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~\cite{MZI_rank}, all target splitting angles $\theta$ are realizable.  For convenience, we derive a quantity called the \textit{coverage}, cov(N)~\cite{self_con_1}, from this probability:
\begin{align}
        \text{cov}(N) &=\prod_n\left(1-P_n(\theta\!<\!\theta_{\rm min})-P_n(\theta\!>\!\theta_{\rm max})\right) \nonumber \\
    	& \approx \exp\Bigl(-\sum_n{\bigl[P_n(\theta\!<\!\theta_{\rm min}) + P_n(\theta\!>\!\theta_{\rm max})\bigr]}\Bigr)
    \label{eq:cov}
\end{align}
For a \texttt{CLEMENTS} mesh, Ref.~\cite{self_con_1} proves that $\text{cov}^{\text{(clem)}}(N) = \text{exp}\left( -\frac{N^{3}\sigma^{2}}{3} \right)$. Like the matrix error in the previous subsection, the coverage of an SCF mesh is computed using Taylor series expansions. The result is:
\begin{equation}
    \text{cov}^{\text{(scf)}}(N) = \exp\Bigl(-\frac{8N^2\log_2(N)}{\pi^2} \sigma^2\Bigr) \label{eq:covft}
\end{equation}
which is greater than $\text{cov}^{\text{(clem)}}(N)$ for the same $\sigma$ for all but very small $N$. 

\section{\label{sec: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.~\cite{NN_arch_1, NN_arch_2, NN_arch_3}\textemdash each neural net layer is implemented by an SCF mesh connected to an electro-optic nonlinearity~\cite{MZI_NL} (Fig.~\ref{fig_5}(a)). All our networks had two layers. Our simulations used the \texttt{meshes}~\cite{meshes_ryan} package, and results are presented for the MNIST~\cite{MNIST} image classification task. The preprocessing of the images involved low-pass filtering and was identical to the procedure adopted in Refs.~\cite{NN_arch_3,hec_1}. The standard Cross Entropy loss and the Adam optimizer were used for training.

ONNs trained with SCF meshes achieved accuracies that matched the \texttt{CLEMENTS} mesh~\cite{hec_1, self_con_1} of $\sim$ 95\%-96\% for small meshes ($N = 64$) and $\sim$ 97\% for larger meshes ($N = 256$) . Next, we simulated the effect of MZI errors on the trained SCF mesh neural net; Fig.~\ref{fig_5}(b) shows the median classification accuracy of 10 independently trained networks as a function of splitter errors. SCF networks of size $N = 64$ yielded $\sim$ 95\% test accuracy while those of size $N = 256$ reached $\sim$ 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 $\sim$ 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~\cite{process_err} and larger than the corresponding cutoff for \texttt{CLEMENTS} meshes (which is 6\%~\cite{hec_1}).

\begin{figure}[t]
    \centering
    \includegraphics[width = \columnwidth]{Figure_pdfs/Fig_4_final.pdf}
    \caption{The SCF mesh at different stages of pruning: \textbf{(a)} A maximally pruned mesh, which is the standard FFT mesh, $\alpha = 1.0$. \textbf{(b)} $\alpha = 1.5$ \textbf{(c)} $\alpha = 1.75$ }
    \label{fig_4}
\end{figure}

\subsection*{Weight pruning}

The number of columns of MZIs with stride $s$ is $\frac{N}{2s}$ 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 \textit{fractal dimension} $\alpha \in [1, 2]$\textemdash in a partially pruned mesh, the number of columns of stride $s$ is $\left( \frac{N}{2s} \right)^{\alpha - 1}$. Setting $\alpha = 1$ and $\alpha = 2$ yields the FFT and SCF meshes respectively. Controlling $\alpha$ allows us to tune the number of degrees of freedom and reduce the areal footprint of the device while ensuring full connectivity. Fig.~\ref{fig_4} illustrates partially pruned $8 \times 8$ SCF meshes for different values of $\alpha$. The depth of the pruned network (approximated to leading order in $N$) scales as $\frac{N^{\alpha - 1}}{2^{\alpha - 1} - 1}$ while the number of MZIs is approximately $\frac{N^{\alpha}}{2\left(2^{\alpha} - 1\right)}$. Care has to be taken during pruning to ensure that no two consecutive columns have the same stride\textemdash any such columns would collapse into a single column.

Fig.~\ref{fig_5}(c) illustrates the results of training networks with pruned meshes of ideal MZIs for different $\alpha$. Increasing $\alpha$ (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 \times 64$ 2-layer ONN already achieves 95\% accuracy, which is commensurate with the performance of present-day DNN accelerators~\cite{Opt_DNN_Ryan}.

\begin{figure}
    \centering
    \includegraphics[width = \columnwidth]{Figure_pdfs/Fig_5_final.pdf}
    \caption{\textbf{(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. \textbf{(b)} Simulated classification accuracy as a function of mesh size ($N =$ 64 and 256) and MZI error $\sigma$. \textbf{(c)} Classification accuracy as a function of degree of pruning $\alpha \in [1, 2]$ and mesh size ($N = $ 16 and 64 and 256). \textbf{(d)}, \textbf{(e)} Classification accuracy as a function of the degree of pruning $\alpha$ and MZI error, trained on maximally faulty mesh of size $N =$ 64 and 256.  }
    \label{fig_5}
\end{figure}

\subsection*{Lower bounds on SCF mesh performance}

We performed error-aware `maximally faulty mesh' training \cite{Sri_paper} to obtain non-trivial (empirical) lower bounds on the performance of SCF mesh ONNs with MZI errors as large as $30\%$. Ref.~\cite{Sri_paper} shows that matrices obtained by training more faulty networks can be exactly transferred to less faulty ones. Maximally faulty 2-layer SCF mesh ONNs of sizes 64 and 256 were trained for several values of MZI error level $\sigma$ and fractal dimension $\alpha$; the resultant test accuracies are depicted in Fig.~\ref{fig_5}(d),(e). As expected, the 256 sized ONNs perform better than the 64 sized ones. However, it is particularly striking in the size 256 case that one can prune the network up to $\alpha\sim1.1$ and still lose only 1\% in classification accuracy. The size 64 case allows less freedom in pruning because it had a lower number of parameters to begin with. The performance also seems to be a nearly constant function of the MZI error. This investigation indicates that one can aggressively prune very faulty ONNs but still achieve excellent performance.   

\section{\label{sec:5}Conclusion}

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)$ vs. $O(\sqrt{N \log_{2}(N)})$ significantly impacts the performance for larger meshes. This reduced scaling with $N$ implies that SCF meshes are more senstive to improvements in foundry process. Reducing $\sigma$ by $2\times$ corresponds to a $4\times$ increase in the maximum mesh size for \textsc{Reck} and \textsc{Clements} owing to the $O(N\sigma^2)$ scaling; for the butterfly fractal, the corresponding increase would be $16\times$.


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:
\begin{itemize}
    \item \textit{Incomplete Decomposition}: To minimize the number of crossings, the block decomposition of the mesh can be terminated such that the smallest block size $N_{\mathrm{blk}} > 2$. Each block will then take the form of a standard \texttt{CLEMENTS} geometry with no intra-block crossings. 
    
    \item \textit{Out-of-plane Crossings}: Inter-block crossings ($s > 1$ among arbitrary block sizes) can be implemented using waveguide escalators and out-of-plane crossings. These crossings have been shown to have much lower crosstalk than in-plane crossings~\cite{3d_cross_1, 3d_cross_2}. 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~\cite{3d_cross_polymer_1, 3d_cross_polymer_2} or glass~\cite{3d_cross_glass_1, 3d_cross_glass_2}, and could utilize ``hockey-stick'' escalator couplers to reduce the alignment tolerances for each chiplet \cite{HockeyStick}. 
\end{itemize}
This suggests that there exists an optimal block or chiplet size $N_{\mathrm{blk}}$ for multi-chip modules, which will have to be determined by the trade-off between intra-block errors (that favors small $N_{\mathrm{blk}}$) and inter-block losses (that favors large $N_{\mathrm{blk}}$).

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 $\mathcal{E}_c^{\mathrm{(clem)}} \gtrsim 0.2$ even in the presence of error correction. This brings us to the regime below the recommended 4--8 bits of precision that is usually targeted for DNN accelerators~\cite{DNN_acc_1, DNN_acc_2, DNN_acc_3}. 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 a frontrunning candidate architecture for future systems.  

\section*{Acknowledgements}

J.R.B.\ would like to acknowledge support from NTT Research, Inc.\ for the duration of this project. R.H., S.B., and D.E.\ would like to thank Prof.\ David A.B.\ Miller and Dr.\ Sunil Pai for helpful discussions. The simulations presented in this paper were performed on the DeepThought2 supercomputing cluster. The authors acknowledge the University of Maryland supercomputing resources (\url{http://hpcc.umd.edu}) made available for conducting the research reported in this paper.

