# A Self-Similar Sine-Cosine Fractal Architecture for Multiport
Interferometers

Jasvith Raj Basani  Department of Electrical and Computer Engineering,
Institute for Research in Electronics and Applied Physics, and Joint Quantum
Institute, University of Maryland, College Park, MD 20742, USA  Sri Krishna
Vadlamani  Research Laboratory of Electronics, Massachusetts Institute of
Technology, 77 Massachusetts Avenue, Cambridge, Massachusetts 02139, USA
Saumil Bandyopadhyay  Research Laboratory of Electronics, Massachusetts
Institute of Technology, 77 Massachusetts Avenue, Cambridge, Massachusetts
02139, USA  Dirk R. Englund  Research Laboratory of Electronics, Massachusetts
Institute of Technology, 77 Massachusetts Avenue, Cambridge, Massachusetts
02139, USA  Ryan Hamerly  Research Laboratory of Electronics, Massachusetts
Institute of Technology, 77 Massachusetts Avenue, Cambridge, Massachusetts
02139, USA  NTT Research Inc., PHI Laboratories, 940 Stewart Drive, Sunnyvale,
CA 94085, USA

(August 27, 2024)

###### Abstract

Multiport interferometers based on integrated beamsplitter meshes have
recently captured interest as a platform for many emerging technologies. In
this paper, we present a novel architecture for multiport interferometers
based on the Sine-Cosine fractal decomposition of a unitary matrix. Our
architecture is unique in that it is self-similar, enabling the construction
of modular multi-chiplet devices. Due to this modularity, our design enjoys
improved resilience to hardware imperfections as compared to conventional
multiport interferometers. Additionally, the structure of our circuit enables
systematic truncation, which is key in reducing the hardware footprint of the
chip as well as compute time in training optical neural networks, while
maintaining full connectivity. Numerical simulations show that truncation of
these meshes gives robust performance even under large fabrication errors.
This design is a step forward in the construction of large-scale programmable
photonics, removing a major hurdle in scaling up to practical machine learning
and quantum computing applications.

††preprint: APS/123-QED

##  I 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, 2, 3, 4, 5], machine learning [6, 7], optimization [8],
and quantum simulation [9, 10, 11, 12, 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-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, 16, 17,
18, 19, 20], local correction [21, 22], self-configuration [23, 24, 25, 26],
and hardware augmentation [27, 28, 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\times 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 II introduces and discusses the
Sine-Cosine Fractal architecture. Section III contains analytical and
numerical results on the expressivity of SCF meshes in the presence of
beamsplitter errors. Section IV reports the performance of ONNs constructed
from both complete and pruned SCF meshes, and Section V concludes the paper
with a further discussion of the scope and impact of our work.

##  II 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. 1(a)) is given by:

$\displaystyle T(\theta,\phi)$ $\displaystyle=\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}$
$\displaystyle=ie^{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}$ (1)

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. 1(a) depicts a mesh that implements an $8\times 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.

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 [30] shows that the block-wise SVD of $U$
yields:

$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}$
(2)

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{bmatrix}D_{11}&D_{12}\\\
D_{21}&D_{22}\end{bmatrix}=ie^{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}$ (3)

with $(\Theta,\Phi)$ representing diagonal matrices that encode phase shifts.
Fig. 1(b) depicts Eq. (2) graphically for the $8\times 8$ case—the $4\times 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}$ subblocks to obtain the mesh of Fig.
1(c); the $2\times 2$ unitaries obtained from the $4\times 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^{n}$) matrix: one
recursively performs blockwise SVDs on each unitary matrix of size greater
than $2$ until the full decomposition consists only of $2\times 2$ matrices
that connect different modes.

Figure 1: (a) $8\times 8$ CLEMENTS mesh. (b) First step of the block
decomposition of the mesh. This results in $4\times 4$ quadrants. (c) Further
decomposition of the quadrants, which results in the SCF mesh.

Like the 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.
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.

##  III 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.

###  III.1 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 [39] induces a distribution $P(\theta)$ over the
phase-shifts. Russell et al. [40] show that the phase-shift distribution of
the $n$-th MZI (according to any indexing scheme) in either the RECK or
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,\ldots,N-1\\}$, called the 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$ [40]. 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. 2).

However, unlike the RECK and 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. (2) 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. 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:

$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}}$
(4)

###  III.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 ($\alpha,\beta$). These splitter errors
perturb the MZI transfer matrix as follows:

$\displaystyle T^{\prime}(\theta^{\prime},\phi^{\prime},\alpha,\beta)$
$\displaystyle=\begin{bmatrix}\cos(\tfrac{\pi}{4}\\!+\\!\beta)&i\sin(\tfrac{\pi}{4}\\!+\\!\beta)\\\
i\sin(\tfrac{\pi}{4}\\!+\\!\beta)&\cos(\tfrac{\pi}{4}\\!+\\!\beta)\end{bmatrix}\begin{bmatrix}e^{i\theta^{\prime}}&0\\\
0&1\end{bmatrix}$ $\displaystyle\ \ \ \
\times\begin{bmatrix}\cos(\tfrac{\pi}{4}\\!+\\!\alpha)&i\sin(\tfrac{\pi}{4}\\!+\\!\alpha)\\\
i\sin(\tfrac{\pi}{4}\\!+\\!\alpha)&\cos(\tfrac{\pi}{4}\\!+\\!\alpha)\end{bmatrix}\begin{bmatrix}e^{i\phi^{\prime}}&0\\\
0&1\end{bmatrix}$ (5)

Bandyopadhyay et al. [21] show that it is always possible to choose phase-
shifts $\theta^{\prime},\phi^{\prime}$ 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:

$\underbrace{2|\alpha+\beta|}_{\theta_{\rm
min}}\leq\theta\leq\underbrace{\pi-2|\alpha-\beta|}_{\theta_{\rm max}}$ (6)

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.

Figure 2: Distribution of crossing angles for the CLEMENTS and SCF meshes as a
function of mesh size and stride respectively. The crossing angles of the
CLEMENTS mesh increasingly clusters near $\theta=0$ for larger meshes.
Crossing angles for the SCF mesh become increasingly uniformly distributed as
the stride increases.

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. [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 $\mathcal{E}_{0}$ (for both CLEMENTS and SCF meshes)
has been shown [24] to scale linearly with the MZI error $\sigma$ as
$\mathcal{E}_{0}=\sqrt{2N}\sigma$. 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, $\mathcal{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):

$\displaystyle(\mathcal{E}_{c})^{2}$
$\displaystyle=\frac{1}{2N}\sum_{n}\Bigl{[}\int_{0}^{\theta_{\rm
min}}P_{n}(\theta)(\theta-\theta_{\rm min})^{2}\mathrm{d}\theta$
$\displaystyle\qquad\qquad\ \ +\int_{\theta_{\rm
max}}^{\pi}P_{n}(\theta)(\theta-\theta_{\rm max})^{2}\mathrm{d}\theta\Bigr{]}$
(7)

Eq. (7) indicates that the effectiveness of error-correction, measured by
$\mathcal{E}_{c}$, is strongly dependent on the distribution $P_{n}(\theta)$.
For the CLEMENTS mesh, Ref. [24] 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. (7) are estimated for the Sine-Cosine Fractal mesh by
Taylor-expanding Eq. (4) to first order about $\theta=0,\pi$ respectively. The
result is:

$\mathcal{E}_{c}^{\mathrm{(scf)}}=\frac{4}{\pi}\sqrt{N\log_{2}(N)}\,\sigma^{2}.$
(8)

The ratio of corrected errors for both meshes is:

$\frac{\mathcal{E}_{c}^{\rm(clem)}}{\mathcal{E}_{c}^{\rm(scf)}}=\sqrt{\frac{\pi^{2}}{24}\frac{N}{\log_{2}(N)}}$
(9)

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

Figure 3: (a) Scaling of matrix error with process variation $\sigma$. (b)
Scaling of matrix error with mesh size, showing the advantage of the SCF mesh
over the CLEMENTS mesh for larger mesh sizes.

We performed numerical experiments on meshes up to size $N=1024$ to validate
the above expressions (Fig. 3). The observed uncorrected error
$\mathcal{E}_{0}$ of both meshes (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 CLEMENTS mesh. 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. 3(b) in the large scale regime ($N>1024$).

###  III.3 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 [40], all target
splitting angles $\theta$ are realizable. For convenience, we derive a
quantity called the coverage, cov(N) [24], from this probability:

$\displaystyle\text{cov}(N)$
$\displaystyle=\prod_{n}\left(1-P_{n}(\theta\\!<\\!\theta_{\rm
min})-P_{n}(\theta\\!>\\!\theta_{\rm max})\right)$
$\displaystyle\approx\exp\Bigl{(}-\sum_{n}{\bigl{[}P_{n}(\theta\\!<\\!\theta_{\rm
min})+P_{n}(\theta\\!>\\!\theta_{\rm max})\bigr{]}}\Bigr{)}$ (10)

For a CLEMENTS mesh, Ref. [24] 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:

$\text{cov}^{\text{(scf)}}(N)=\exp\Bigl{(}-\frac{8N^{2}\log_{2}(N)}{\pi^{2}}\sigma^{2}\Bigr{)}$
(11)

which is greater than $\text{cov}^{\text{(clem)}}(N)$ for the same $\sigma$
for all but very small $N$.

##  IV 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 [41] (Fig. 5(a)). All our
networks had two layers. Our simulations used the meshes [42] package, and
results are presented for the MNIST [43] image classification task. The
preprocessing of the images involved low-pass filtering and was identical to
the procedure adopted in Refs. [26, 21]. 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 $\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. 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 [44] and larger than the
corresponding cutoff for CLEMENTS meshes (which is 6% [21]).

Figure 4: The SCF mesh at different stages of pruning: (a) A maximally pruned
mesh, which is the standard FFT mesh, $\alpha=1.0$. (b) $\alpha=1.5$ (c)
$\alpha=1.75$

### 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 fractal dimension $\alpha\in[1,2]$—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. 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—any such columns would collapse into a single column.

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 [45].

Figure 5: (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 $\sigma$. (c)
Classification accuracy as a function of degree of pruning $\alpha\in[1,2]$
and mesh size ($N=$ 16 and 64 and 256). (d), (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.

### Lower bounds on SCF mesh performance

We performed error-aware ‘maximally faulty mesh’ training [46] to obtain non-
trivial (empirical) lower bounds on the performance of SCF mesh ONNs with MZI
errors as large as $30\%$. Ref. [46] 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. 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\sim
1.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.

##  V 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 Reck and
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:

  * •

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
CLEMENTS geometry with no intra-block crossings.

  * •

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 [47, 48]. 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 [49, 50] or glass [51, 52], and could utilize
“hockey-stick” escalator couplers to reduce the alignment tolerances for each
chiplet [53].

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 [54, 55, 56]. 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.

## 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 (http://hpcc.umd.edu) made available for conducting the research
reported in this paper.

## References

  * Annoni  _et al._ [2017] A. Annoni, E. Guglielmi, M. Carminati, G. Ferrari, M. Sampietro, D. A. Miller, A. Melloni, and F. Morichetti, Unscrambling light—automatically undoing strong mixing between modes, Light: Science & Applications 6, e17110 (2017).
  * Ribeiro  _et al._ [2016] A. Ribeiro, A. Ruocco, L. Vanacker, and W. Bogaerts, Demonstration of a 4$\times$ 4-port universal linear circuit, Optica 3, 1348 (2016).
  * Milanizadeh  _et al._ [2019] M. Milanizadeh, P. Borga, F. Morichetti, D. Miller, and A. Melloni, Manipulating free-space optical beams with a silicon photonic mesh, in  _2019 IEEE Photonics Society Summer Topical Meeting Series (SUM)_ (IEEE, 2019) pp. 1–2.
  * Zhuang  _et al._ [2015] L. Zhuang, C. G. Roeloffzen, M. Hoekman, K.-J. Boller, and A. J. Lowery, Programmable photonic signal processor chip for radiofrequency applications, Optica 2, 854 (2015).
  * Notaros  _et al._ [2017] J. Notaros, J. Mower, M. Heuck, C. Lupo, N. C. Harris, G. R. Steinbrecher, D. Bunandar, T. Baehr-Jones, M. Hochberg, S. Lloyd, _et al._ , Programmable dispersion on a photonic integrated circuit for classical and quantum applications, Optics Express 25, 21275 (2017).
  * Shen  _et al._ [2017] Y. Shen, N. C. Harris, S. Skirlo, M. Prabhu, T. Baehr-Jones, M. Hochberg, X. Sun, S. Zhao, H. Larochelle, D. Englund, _et al._ , Deep learning with coherent nanophotonic circuits, Nature Photonics 11, 441 (2017).
  * Basani  _et al._ [2022] J. R. Basani, M. Heuck, D. R. Englund, and S. Krastanov, All-photonic artificial neural network processor via non-linear optics, arXiv preprint arXiv:2205.08608 arXiv:2205.08608 (2022).
  * Prabhu  _et al._ [2020] M. Prabhu, C. Roques-Carmes, Y. Shen, N. Harris, L. Jing, J. Carolan, R. Hamerly, T. Baehr-Jones, M. Hochberg, V. Čeperić, _et al._ , Accelerating recurrent ising machines in photonic integrated circuits, Optica 7, 551 (2020).
  * Harris  _et al._ [2017] N. C. Harris, G. R. Steinbrecher, M. Prabhu, Y. Lahini, J. Mower, D. Bunandar, C. Chen, F. N. Wong, T. Baehr-Jones, M. Hochberg, _et al._ , Quantum transport simulations in a programmable nanophotonic processor, Nature Photonics 11, 447 (2017).
  * Wang  _et al._ [2018] J. Wang, S. Paesani, Y. Ding, R. Santagati, P. Skrzypczyk, A. Salavrakos, J. Tura, R. Augusiak, L. Mančinska, D. Bacco, _et al._ , Multidimensional quantum entanglement with large-scale integrated optics, Science 360, 285 (2018).
  * Qiang  _et al._ [2018] X. Qiang, X. Zhou, J. Wang, C. M. Wilkes, T. Loke, S. O’Gara, L. Kling, G. D. Marshall, R. Santagati, T. C. Ralph, _et al._ , Large-scale silicon quantum photonics implementing arbitrary two-qubit processing, Nature Photonics 12, 534 (2018).
  * Sparrow  _et al._ [2018] C. Sparrow, E. Martín-López, N. Maraviglia, A. Neville, C. Harrold, J. Carolan, Y. N. Joglekar, T. Hashimoto, N. Matsuda, J. L. O’Brien, _et al._ , Simulating the vibrational quantum dynamics of molecules using photonics, Nature 557, 660 (2018).
  * Carolan  _et al._ [2015] J. Carolan, C. Harrold, C. Sparrow, E. Martín-López, N. J. Russell, J. W. Silverstone, P. J. Shadbolt, N. Matsuda, M. Oguma, M. Itoh, _et al._ , Universal linear optics, Science 349, 711 (2015).
  * Fang  _et al._ [2019] M. Y.-S. Fang, S. Manipatruni, C. Wierzynski, A. Khosrowshahi, and M. R. DeWeese, Design of optical neural networks with component imprecisions, Optics Express 27, 14009 (2019).
  * Burgwal  _et al._ [2017] R. Burgwal, W. R. Clements, D. H. Smith, J. C. Gates, W. S. Kolthammer, J. J. Renema, and I. A. Walmsley, Using an imperfect photonic network to implement random unitaries, Optics Express 25, 28236 (2017).
  * Mower  _et al._ [2015] J. Mower, N. C. Harris, G. R. Steinbrecher, Y. Lahini, and D. Englund, High-fidelity quantum state evolution in imperfect photonic integrated circuits, Physical Review A 92, 032322 (2015).
  * López [2019] D. P. López, Programmable integrated silicon photonics waveguide meshes: optimized designs and control algorithms, IEEE Journal of Selected Topics in Quantum Electronics 26, 1 (2019).
  * López  _et al._ [2020] A. López, D. Pérez, P. DasMahapatra, and J. Capmany, Auto-routing algorithm for field-programmable photonic gate arrays, Optics Express 28, 737 (2020).
  * Pérez-López  _et al._ [2020] D. Pérez-López, A. López, P. DasMahapatra, and J. Capmany, Multipurpose self-configuration of programmable photonic circuits, Nature Communications 11, 1 (2020).
  * Pai  _et al._ [2019] S. Pai, B. Bartlett, O. Solgaard, and D. A. Miller, Matrix optimization on universal unitary photonic devices, Physical Review Applied 11, 064044 (2019).
  * Bandyopadhyay  _et al._ [2021] S. Bandyopadhyay, R. Hamerly, and D. Englund, Hardware error correction for programmable photonics, Optica 8, 1247 (2021).
  * Kumar  _et al._ [2021] S. P. Kumar, L. Neuhaus, L. G. Helt, H. Qi, B. Morrison, D. H. Mahler, and I. Dhand, Mitigating linear optics imperfections via port allocation and compilation, arXiv preprint arXiv:2103.03183 arXiv:2103.03183 (2021).
  * Miller [2017] D. A. Miller, Setting up meshes of interferometers–reversed local light interference method, Optics Express 25, 29233 (2017).
  * Hamerly  _et al._ [2022a] R. Hamerly, S. Bandyopadhyay, and D. Englund, Stability of self-configuring large multiport interferometers, Physical Review Applied 18, 024018 (2022a).
  * Hamerly  _et al._ [2022b] R. Hamerly, S. Bandyopadhyay, and D. Englund, Accurate self-configuration of rectangular multiport interferometers, Physical Review Applied 18, 024019 (2022b).
  * Pai  _et al._ [2020] S. Pai, I. A. Williamson, T. W. Hughes, M. Minkov, O. Solgaard, S. Fan, and D. A. Miller, Parallel programming of an arbitrary feedforward photonic network, IEEE Journal of Selected Topics in Quantum Electronics 26, 1 (2020).
  * Hamerly  _et al._ [2021] R. Hamerly, S. Bandyopadhyay, and D. Englund, Infinitely scalable multiport interferometers, arXiv preprint arXiv:2109.05367 arXiv:2109.05367 (2021).
  * Suzuki  _et al._ [2015] K. Suzuki, G. Cong, K. Tanizawa, S.-H. Kim, K. Ikeda, S. Namiki, and H. Kawashima, Ultra-high-extinction-ratio 2$\times$ 2 silicon optical switch with variable splitter, Optics Express 23, 9086 (2015).
  * Miller [2015] D. A. Miller, Perfect optics with imperfect components, Optica 2, 747 (2015).
  * Polcari [2018] J. Polcari, Generalizing the Butterfly Structure of the FFT, in  _Advanced Research in Naval Engineering_ (Springer, 2018) pp. 35–52.
  * Reck  _et al._ [1994] M. Reck, A. Zeilinger, H. J. Bernstein, and P. Bertani, Experimental realization of any discrete unitary operator, Physical Review Letters 73, 58 (1994).
  * Clements  _et al._ [2016] W. R. Clements, P. C. Humphreys, B. J. Metcalf, W. S. Kolthammer, and I. A. Walmsley, Optimal design for universal multiport interferometers, Optica 3, 1460 (2016).
  * Cooley and Tukey [1965] J. W. Cooley and J. W. Tukey, An algorithm for the machine calculation of complex Fourier series, Mathematics of computation 19, 297 (1965).
  * Flamini  _et al._ [2017] F. Flamini, N. Spagnolo, N. Viggianiello, A. Crespi, R. Osellame, and F. Sciarrino, Benchmarking integrated linear-optical architectures for quantum information processing, Scientific Reports 7, 1 (2017).
  * Gu  _et al._ [2020] J. Gu, Z. Zhao, C. Feng, M. Liu, R. T. Chen, and D. Z. Pan, Towards area-efficient optical neural networks: an FFT-based architecture, in  _2020 25th Asia and South Pacific Design Automation Conference (ASP-DAC)_ (IEEE, 2020) pp. 476–481.
  * Hinton  _et al._ [2012] G. E. Hinton, N. Srivastava, A. Krizhevsky, I. Sutskever, and R. R. Salakhutdinov, Improving neural networks by preventing co-adaptation of feature detectors, arXiv preprint arXiv:1207.0580 arXiv:1207.0580 (2012).
  * Tanaka  _et al._ [2020] H. Tanaka, D. Kunin, D. L. Yamins, and S. Ganguli, Pruning neural networks without any data by iteratively conserving synaptic flow, Advances in Neural Information Processing Systems 33, 6377 (2020).
  * Lee and Park [2007] H.-Y. Lee and I.-C. Park, Balanced binary-tree decomposition for area-efficient pipelined FFT processing, IEEE Transactions on Circuits and Systems I: Regular Papers 54, 889 (2007).
  * Tung [2003] W.-K. Tung, Group Theory in Physics: An Introduction to Symmetry Principles, Group Representations, and Special Functions in Classical and Quantum Physics (2003).
  * Russell  _et al._ [2017] N. J. Russell, L. Chakhmakhchyan, J. L. O’Brien, and A. Laing, Direct dialling of Haar random unitary matrices, New Journal of Physics 19, 033007 (2017).
  * Williamson  _et al._ [2019] I. A. Williamson, T. W. Hughes, M. Minkov, B. Bartlett, S. Pai, and S. Fan, Reprogrammable electro-optic nonlinear activation functions for optical neural networks, IEEE Journal of Selected Topics in Quantum Electronics 26, 1 (2019).
  * Hamerly [2021] R. Hamerly, Meshes: Tools for modeling photonic beamsplitter mesh networks, https://github.com/QPG-MIT/meshes (2021).
  * LeCun [1998] Y. LeCun, The MNIST database of handwritten digits, http://yann.lecun.com/exdb/mnist/ (1998).
  * Mikkelsen  _et al._ [2014] J. C. Mikkelsen, W. D. Sacher, and J. K. Poon, Dimensional variation tolerant silicon-on-insulator directional couplers, Optics Express 22, 3145 (2014).
  * Hamerly  _et al._ [2019] R. Hamerly, L. Bernstein, A. Sludds, M. Soljačić, and D. Englund, Large-scale optical neural networks based on photoelectric multiplication, Physical Review X 9, 021032 (2019).
  * Vadlamani  _et al._ [2022] S. K. Vadlamani, D. Englund, and R. Hamerly, Transferable Learning on Analog Hardware, In preparation (2022).
  * Jones  _et al._ [2013] A. M. Jones, C. T. DeRose, A. L. Lentine, D. C. Trotter, A. L. Starbuck, and R. A. Norwood, Ultra-low crosstalk, CMOS compatible waveguide crossings for densely integrated photonic interconnection networks, Optics Express 21, 12002 (2013).
  * Sacher  _et al._ [2018] W. D. Sacher, J. C. Mikkelsen, Y. Huang, J. C. Mak, Z. Yong, X. Luo, Y. Li, P. Dumais, J. Jiang, D. Goodwill, _et al._ , Monolithically integrated multilayer silicon nitride-on-silicon waveguide platforms for 3-D photonic circuits and devices, Proceedings of the IEEE 106, 2232 (2018).
  * Lindenmann  _et al._ [2012] N. Lindenmann, G. Balthasar, D. Hillerkuss, R. Schmogrow, M. Jordan, J. Leuthold, W. Freude, and C. Koos, Photonic wire bonding: a novel concept for chip-scale interconnects, Optics Express 20, 17667 (2012).
  * Billah  _et al._ [2018] M. R. Billah, M. Blaicher, T. Hoose, P.-I. Dietrich, P. Marin-Palomo, N. Lindenmann, A. Nesic, A. Hofmann, U. Troppenz, M. Moehrle, _et al._ , Hybrid integration of silicon photonics circuits and InP lasers by photonic wire bonding, Optica 5, 876 (2018).
  * Mennea  _et al._ [2018] P. L. Mennea, W. R. Clements, D. H. Smith, J. C. Gates, B. J. Metcalf, R. H. Bannerman, R. Burgwal, J. J. Renema, W. S. Kolthammer, I. A. Walmsley, _et al._ , Modular linear optical circuits, Optica 5, 1087 (2018).
  * Szameit and Nolte [2010] A. Szameit and S. Nolte, Discrete optics in femtosecond-laser-written photonic structures, Journal of Physics B: Atomic, Molecular and Optical Physics 43, 163001 (2010).
  * Bandyopadhyay and Englund [2021] S. Bandyopadhyay and D. Englund, Alignment-free photonic interconnects, arXiv preprint arXiv:2110.12851 arXiv:2110.12851 (2021).
  * Friedmann  _et al._ [2013] S. Friedmann, N. Frémaux, J. Schemmel, W. Gerstner, and K. Meier, Reward-based learning under hardware constraints—using a RISC processor embedded in a neuromorphic substrate, Frontiers in Neuroscience 7, 160 (2013).
  * Akopyan  _et al._ [2015] F. Akopyan, J. Sawada, A. Cassidy, R. Alvarez-Icaza, J. Arthur, P. Merolla, N. Imam, Y. Nakamura, P. Datta, G.-J. Nam, _et al._ , Truenorth: Design and tool flow of a 65 mw 1 million neuron programmable neurosynaptic chip, IEEE transactions on computer-aided design of integrated circuits and systems 34, 1537 (2015).
  * Jouppi  _et al._ [2017] N. P. Jouppi, C. Young, N. Patil, D. Patterson, G. Agrawal, R. Bajwa, S. Bates, S. Bhatia, N. Boden, A. Borchers, _et al._ , In-datacenter performance analysis of a tensor processing unit, in  _Proceedings of the 44th annual international symposium on computer architecture_ (2017) pp. 1–12.


11institutetext: College of Computing and Informatics, Drexel University,
Philadelphia, PA, USA
11email<EMAIL_ADDRESS>
# Exploring Wasserstein Distance across Concept Embeddings for Ontology
Matching

Yuan An  Alex Kalinowski  Jane Greenberg

###### Abstract

Measuring the distance between ontological elements is fundamental for
ontology matching. String-based distance metrics are notorious for shallow
syntactic matching. In this exploratory study, we investigate Wasserstein
distance targeting continuous space that can incorporate various types of
information. We use a pre-trained word embeddings system to embed ontology
element labels. We examine the effectiveness of Wasserstein distance for
measuring similarity between ontologies, and discovering and refining
matchings between individual elements. Our experiments with the OAEI
conference track and MSE benchmarks achieved competitive results compared to
the leading systems.

###### Keywords:

ontology matching optimal transport Wasserstein Distance ontology embedding.

##  1 Introduction

Semantic and structural heterogeneity is widespread among ontologies. To
bridge the heterogeneity, almost all ontology matching systems [ontology-
matching-literature-review, AML-ontology-matching, LogMap] need to evaluate
the distances between ontological elements. String-based distance metrics have
dominated the field [string-similarity-ontology-matching]. However, string-
based distance metrics are notorious for shallow syntactic matching. It is
also challenging to determine which string similarity measures to use and how
to effectively combine them in matching systems [string-similarity-ontology-
matching]. Aligning elements across different ontologies essentially involves
measuring the semantic similarity/distance of the ontological elements
representing items in the underlying domains. Recently, word embeddings
[word2vec] have been used to successfully encode syntactic and semantic word
relationships. As a result, word embeddings have displayed excellent
performance in applications for cross-lingual word alignment [unsupervised-
alignment-wasserstein]. Ontology matching techniques have also been created
using embeddings, with embedding vectors predominantly utilized as inputs in
supervised [ontology-alignment-embedding-random-forest, BERTMap] or distantly
supervised [augmenting-ontology-alignment] machine learning models.

While machine learning is effective in making use of ontology embeddings, a
significant effort is required to gather training instances. Although some
systems directly use the cosine similarity between embedding vectors
[biomedical-ontology-alignment] for deriving candidate matchings, the
embeddings are first retrofitted through tailored training instances. In
practice, the present top matching algorithms are largely unsupervised and
rely only upon existing ontology and external sources. Motivated by the
desired property of unsupervised learning, we posit the following question:
_If ontology elements can be readily encoded as embedding vectors, can we use
them in an unsupervised fashion for ontology matching?_

To research the question, we formulate ontology matching as an optimal
transport problem from a source ontology embedding space to a target ontology
embedding space. Optimal Transport (OT) [ot-vallani] has been applied to
various alignment applications including word embeddings alignment [Gromov-W],
sequence-to-sequence learning [chen2019improving], heterogeneous domain
alignment [got], and graph comparison and matching [an-optimal-graph-
comparison]. A desired advantage of OT-based approach is _unsupervised
learning_. The OT solution establishes optimal mappings and a shape-based
distance called _Wasserstein distance_ between distributions. In this study,
we explore the effectiveness of Wasserstein distance for ontology matching.
Our inquiry focuses on answering the following 3 questions:

  1. 1.

_How effective is Wasserstein distance for measuring similarity between
(blocks of) ontologies?_

  2. 2.

_How effetive is the coupling matrix accompanying a Wasserstain distance for
deriving alignments between individual ontological elements?_

  3. 3.

_How effective is Wasserstain distance for refining matching candidates?_

The rest of the paper presents our research and is organized as follows.
Section 2 introduces optimal transport and Wasswerstein distance. Section 3
describes the ontology embeddings used in this study. Section 4 measures
ontology similarity in Wasswerstein distance. Section 5 derives matching
candidates. Section LABEL:sec:refining-matching refines matching candidates.
Section LABEL:sec:experimental-results presents our experiment and results.
Section LABEL:sec:discussion comments on the results. Section
LABEL:sec:related-work discusses related work, and finally, Section
LABEL:sec:conclusion concludes the paper with future directions.

##  2 Optimal Transport and Wasserstein Distance

Optimal transport (OT) [ot-vallani] originated as a solution to the problem of
transporting masses from one configuration onto another with _the least
effort_. OT had deep connections with major economic problems, for example, in
logistics, production planning, and network routing, etc. Since then, optimal
transport has been generalized from practical concerns to powerful
mathematical tools for comparing distributions. In particular, given two point
sets modeled as two discrete distributions, optimal transport is an effective
approach for discovering a minimum cost mapping between the two sets.
Considering the set of embedding vectors of an ontology as a set of points, we
can apply optimal transport to discover a minimum cost mapping between two
sets of ontology embeddings.

Formally, given two sets of ontology embeddings
$\mathbf{X}=\\{\mathbf{x}_{i}\in\mathbb{R}^{d},i=1..n\\}$ and
$\mathbf{Y}=\\{\mathbf{y}_{j}\in\mathbb{R}^{d},j=1..m\\}$, where each
embedding is represented as a vector $\mathbf{x}_{i}$ or
$\mathbf{y}_{j}\in\mathbb{R}^{d}$. Let
$\mu=\sum_{i=1}^{n}{p}(\mathbf{x}_{i})\delta_{\mathbf{x}_{i}}$ be the
probability distribution defined on the set $\mathbf{X}$, where
$\delta_{\mathbf{x}_{i}}$ is the Dirac at the point $\mathbf{x}_{i}$ and
$p(\mathbf{x}_{i})$ is a probability weight associated with the point
$\mathbf{x}_{i}$. Similarly, let
$\nu=\sum_{j=1}^{n}{q}(\mathbf{y}_{j})\delta_{\mathbf{y}_{j}}$ be the
probability distribution defined on the set $\mathbf{Y}$, where
$\delta_{\mathbf{y}_{j}}$ is the Dirac at the point $\mathbf{y}_{j}$ and
$q(\mathbf{y}_{j})$ is a probability weight associated with the point
$\mathbf{y}_{j}$. Usually, we consider uniform weights, e.g.,
$p(\mathbf{x}_{i})=\frac{1}{n}$, for $i=1..n$, and
$q(\mathbf{y}_{j})=\frac{1}{m}$, for $j=1..m$. However, if additional
information is provided, $p(\mathbf{x}_{i})$ and $q(\mathbf{y}_{j})$ can
incorporate the information as non-uniform distributions. Optimal transport
(OT) defines an optimal plan for mass transportation and a distance between
the two distributions.

Specifically, let $\mathbf{C}=[c(\mathbf{x}_{i},\mathbf{y}_{j})]_{i,j}$ be a
ground cost matrix with $c(\mathbf{x}_{i},\mathbf{y}_{j})$ measuring a ground
distance between the individual embeddings $\mathbf{x}_{i}$ and
$\mathbf{y}_{j}$. Let $\mathbf{T}=[T(\mathbf{x}_{i},\mathbf{y}_{j})]_{i,j}$ be
a matrix of a transport plan (or _couplings_) with
$T(\mathbf{x}_{i},\mathbf{y}_{j})$ specifying how much mass will be
transported from point $\mathbf{x}_{i}$ to point $\mathbf{y}_{j}$. Let
$\Pi(\mu,\nu)$ be the set of all feasible transport plans defined as:

$\Pi(\mu,\nu)\stackrel{{\scriptstyle\text{def}}}{{=}}\\{\mathbf{T}\in\mathbb{R}_{+}^{n\times
m}|\mathbf{T}\mathds{1}_{n}=\mu,\mathbf{T}^{\top}\mathds{1}_{m}=\nu\\}$ (1)

where $\mathds{1}_{n}$ and $\mathds{1}_{m}$ are all one vectors,
$\mathbf{T}\mathds{1}_{n}=\mu$ and $\mathbf{T}^{\top}\mathds{1}_{m}=\nu$ are
marginal constraints on feasible plans. The _Optimal Transport_ problem is to
find the map $\mathbf{T}:\mathbf{X}\rightarrow\mathbf{Y}$, where

$\mathbf{T}=\operatorname*{argmin}_{\mathbf{T}\in\Pi(\mu,\nu)}\sum_{i=1}^{n}\sum_{j=1}^{m}c(\mathbf{x}_{i},\mathbf{y}_{j})\cdot
T(\mathbf{x}_{i},\mathbf{y}_{j}),\text{s.t.},\mathbf{T}\mathds{1}_{n}=\mu,\mathbf{T}^{\top}\mathds{1}_{m}=\nu$
(2)

The map $\mathbf{T}$ is also called _coupling matrix_. It gives rise to a
distance measure between the two distributions called _Wasserstein distance_
defined as:

$W(\mu,\nu)\stackrel{{\scriptstyle\text{def}}}{{=}}\min_{\mathbf{T}\in\Pi(\mu,\nu)}\langle\mathbf{C},\mathbf{T}\rangle{=}\min_{\mathbf{T}\in\Pi(\mu,\nu)}\sum_{i=1}^{n}\sum_{j=1}^{m}c(\mathbf{x}_{i},\mathbf{y}_{j})\cdot
T(\mathbf{x}_{i},\mathbf{y}_{j})$ (3)

The optimization problem can be efficiently solved by replacing the objective
with an entropy regularized objective such as in the sinkhorn algorithm [ot-
vallani].

Figure 1: Optimal Transport Problem between a Source and Target Ontology

###### Example 1

_Figure 1 illustrates the configuration of the optimal transport problem
between a source and target ontology. The source ontology has 3 concepts:
Document, Writer and Topic. The target ontology has 4 concepts which are:
Paper, Accepted_Paper, Author and Subject_Area. The table in Figure 1 (a)
shows the Euclidean distances as the ground costs between the embeddings of
the source and target ontology concepts. The points in both the source and
target embedding sets are uniformly distributed, as reported in the last
column $\mu$ and the last row $\nu$. By solving Equation (2), the table in
Figure 1 (b) displays the optimal transport couplings between the points in
the source and target embedding spaces. Notice the couplings are feasible
because they satisfy the following marginal constraints:
$\mathbf{T}\mathds{1}_{n}=\mu$ and $\mathbf{T}^{\top}\mathds{1}_{m}=\nu$. By
solving the Equation (3), the Wasserstein distance between the two embedding
sets is 0.47 given the ground costs and the couplings. _

##  3 Ontology Embeddings and Matching

Ontology embedding is the problem of encoding ontological elements as
numerical vectors. Various methods have been proposed for representing
individual components in an ontology as embeddings. For example,
translational-based methods [transe, transr] and graph neural networks (GNN)
[graph-representation-learning] encode an ontology based on its graph
structure. Text-enhanced methods for ontology embeddings [utilizing-textual,
unsupervised-embedding-enhancements, an-etal-2018-accurate, model-text-
enhanced, Wang2016TextEnhancedRL] encode lexical words of ontology elements.
Logic-aware methods [rocktaschel-etal-2015-injecting, demeester-
etal-2016-lifted] incorporate logical constraints into ontology embeddings. It
is attempting to encode ontologies using the above methods and directly apply
optimal transport for discovering mapping. However, simply applying these
embeddings for ontology matching is negatively impacted by the _lack of
registration_ problem [pmlr-v108-alvarez-melis20a]. Specifically, these
embeddings were mainly developed for applications concerned with a single
ontology, for example, link prediction [open-challenge-inductive-link-
prediction]. The embedding spaces of two independent ontologies may mismatch
due to different dimensions or various rotations and translations. As a
result, there may not exist a distance or the direct geometric locations
between the points in the embedding spaces may not reflect their underlying
genuine relationships. This is a significant issue for optimal transport-based
approach which needs a meaningful ground cost.

We will study the problem of matching the embeddings incorporating structural
and logical information in future work. In this exploratory study, we focus on
the ontology embeddings corresponding to the labels of the ontology elements,
for example, the _‘rdfs:label’_ of a concept. we apply the pre-trained
language model, _fasttext_ [fasttext2016enriching], to encode the labels of
the set of ontology concepts, object and datatype properties. Using the same
pre-trained language model to encode the labels of different ontologies will
alleviate the lack of registration problem, because the resultant embeddings
are in the same embedding space. We will show that OT on label embeddings
already produce promising results. For each element, we first normalizes the
element’s label via a sequence of standard text processing steps. If
necessary, the labels are augmented with synonyms. We then split the
normalized label into individual words which in turn are fed into the pre-
trained language model to obtain their corresponding word embeddings. For the
element, we obtain its embedding by computing the average of the set of
embeddings or use the entire set of the embeddings of the individual words for
next-step processing.

##  4 Wasserstain Distance for Measuring Ontology Similarity

Prior to applying optimal transport and Wasserstein distance for ontology
matching, we first analyze some properties of Wasserstein distance for
capturing ontology similarity. A desirable property is that Wasserstein
distance should closely correlate with ontology similarity. That is, the more
similar two ontologies are, the shorter the Wasserstein distance between them
is. Quantitatively measuring the similarity between two ontologies is a very
challenging problem. One option is to count the minimum number of graph edit
operations to transform one graph to another. It is known that graph edit
distance is NP-hard and it depends on a set of graph edit operations. Another
option is the Jaccard Index on sets. If matchings between ontologies are
available, we can define the Jaccard similarity between ontologies as follows:

$Jac(\mathcal{O}_{S},\mathcal{O}_{T})=\frac{|M|}{|\mathcal{O}_{S}|+|\mathcal{O}_{T}|-|M|},$

where $\mathcal{O}_{S}$ and $\mathcal{O}_{T}$ are source and target
ontologies, $|\mathcal{O}_{S}|$ and $|\mathcal{O}_{T}|$ return the number of
concepts in each ontology, $M$ is the set of matchings between
$\mathcal{O}_{S}$ and $\mathcal{O}_{T}$, and $|M|$ gives the number of
matchings.

Figure 2: Regression plot between Wasserstein similarity ($e^{-wd}$) and
Jaccard similarity based on the given matching references for the conference
ontologies

The OAEI Campaign provides a list of ontology matching tasks each with a set
of curated matching references. In particular, the Conference track contains
21 matching cases among 7 conference ontologies. For each case, we compute:
(1) its Jaccard similarity based on the given matching references, and (2) the
Wasserstein distance between the ontology embeddings. We convert a Wasserstein
distance, $wd$, to a Wasserstein similarity, $ws$ as $ws=e^{-wd}$. Figure 2
shows the regression plot between the Wasserstein similarities and Jaccard
similarities. Furthermore, we calculated the Pearson correlation coefficient
(PCC) between the two sets of similarities defined as
$\rho=\frac{cov(X,Y)}{\sigma_{X}\sigma_{Y}}$. The PCC is 0.77, indicating
Wasserstein distance is highly correlated with the Jaccard coefficients when
matchings between two ontologies are known. As a result, Wasserstein distance
exhibits the desirable property for capturing ontology similarity. We can then
leverage this good property for developing algorithms that derive and refine
matchings using optimal transport as presented below.

##  5 Deriving Matching Candidates from Global Coupling Matrix

Given a source $\mathcal{O}_{S}$ and a target $\mathcal{O}_{T}$ ontology, we
first apply optimal transport _globally_ to the ontologies for generating a
set of candidate matchings. We then compute contextual _local_ Wasserstein
distances to refine the matchings by filtering out false positives.

In this section, we describe the building blocks for deriving matching
candidates globally. Let
$\mathbf{X}=\\{\mathbf{x}_{i}\in\mathbb{R}^{d},i=1..n\\}$ be the source
ontology concept embeddings and
$\mathbf{Y}=\\{\mathbf{y}_{j}\in\mathbb{R}^{d},j=1..m\\}$ be the target
ontology concept embeddings. The optimal transport problem defined in Equation
(2) requires as input ground costs
$\mathbf{C}=[c(\mathbf{x}_{i},\mathbf{y}_{j})]_{i,j}$ and probability
distributions $p(\mathbf{x}_{i})$ and $q(\mathbf{y}_{j})$. For label embedding
spaces, we use the Euclidean distances as the ground costs. For the
probability distributions, we estimate non-uniform weights using the shortest
distances between source and target concept embeddings. In particular, for a
source point $\mathbf{x}_{i}\in\mathbf{X}$, let
$d_{i}=\min_{j=1}^{m}c(\mathbf{x}_{i},\mathbf{y}_{j})$ be the shortest
distance from $\mathbf{x}_{i}$ to all embedding points in the target space.
The distribution $p(\mathbf{x}_{i})$ will be inversely proportional to $d_{i}$
for $i=1..n$. In other words, the greater the shortest distance from a source
point to all target points, the less the weight associated with the source
concept. Similarly, we estimate non-uniform probability distribution
$q(\mathbf{y}_{j})$ associated with the target concepts using the shortest
distances from target points to source points.

The solution, $\mathbf{T}=[T(\mathbf{x}_{i},\mathbf{y}_{j})]_{i,j}$,
$i=1..n,j=1..m$, is a coupling matrix between every source and target
embedding point. We test the following two methods for deriving candidate
matchings from the coupling matrix:

  * •

Mutual Nearest Neighbor (MNN): for a
$\mathbf{x}_{p}\in$$\\{\mathbf{x}_{i}\in\mathbb{R}^{d},i=1..n\\}$, find
$\mathbf{y}_{q}\in$$\\{\mathbf{y}_{j}\in\mathbb{R}^{d},j=1..m\\}$, such that,
$T(\mathbf{x}_{p},\mathbf{y}_{q})$ $=$
$\max\\{T(\mathbf{x}_{p},\mathbf{y}_{j}),j=1..m\\}$ and
$T(\mathbf{x}_{p},\mathbf{y}_{q})$ $=$
$\max\\{T(\mathbf{x}_{i},\mathbf{y}_{q}),i=1..n\\}$.

  * •

Top-K Targets (TopK): for a
$\mathbf{x}_{p}\in$$\\{\mathbf{x}_{i}\in\mathbb{R}^{d},i=1..n\\}$, find $k$
targets
$\\{\mathbf{y}_{q_{1}},\mathbf{y}_{q_{2}},..,\mathbf{y}_{q_{k}}\\}\subset$
$\\{\mathbf{y}_{j}\in\mathbb{R}^{d},j=1..m\\}$, such that,
$T(\mathbf{x}_{p},\mathbf{y}_{q_{z}})$ $\geq$
$\max\\{T(\mathbf{x}_{p},\mathbf{y}_{j}),j\neq q_{1}..q_{k}\\}$, for $z=1..k$.