Interferometric Neural Networks
On the one hand, artificial neural networks have many successful applications
in the field of machine learning and optimization. On the other hand,
interferometers are integral parts of any field that deals with waves such as
optics, astronomy, and quantum physics. Here, we introduce neural networks
composed of interferometers and then build generative adversarial networks from
them. Our networks do not have any classical layer and can be realized on
quantum computers or photonic chips. We demonstrate their applicability for
combinatorial optimization, image classification, and image generation. For
combinatorial optimization, our network consistently converges to the global
optimum or remains within a narrow range of it. In multi-class image
classification tasks, our networks achieve accuracies of 93% and 83%. Lastly,
we show their capability to generate images of digits from 0 to 9 as well as
human faces.
Arun Sehrawat
# Interferometric Neural Networks

###### Abstract

On the one hand, artificial neural networks have many successful applications in the field of machine learning and optimization. On the other hand, interferometers are integral parts of any field that deals with waves such as optics, astronomy, and quantum physics. Here, we introduce neural networks composed of interferometers and then build generative adversarial networks from them. Our networks do not have any classical layer and can be realized on quantum computers or photonic chips. We demonstrate their applicability for combinatorial optimization, image classification, and image generation. For combinatorial optimization, our network consistently converges to the global optimum or remains within a narrow range of it. In multi-class image classification tasks, our networks achieve accuracies of 93% and 83%. Lastly, we show their capability to generate images of digits from 0 to 9 as well as human faces.

## I Introduction

Artificial neural networks (NNs) [1; 2; 3] have demonstrated remarkable success in a wide range of applications in various domains such as autonomous vehicles, healthcare, finance, robotics, gaming, over-the-top media services, and chatbots. Some of the applications include image classification [4; 5; 6], natural language processing, speech recognition, robotic control systems, recommendation systems, and image generation [7; 8; 9; 10]. Convolutional NNs (CNNs) were introduced in [4; 5] and further improved in seminal works like [6] for image classification. Similarly, generative adversarial networks (GANs) made their debut in [7] and underwent notable advancements in subsequent works, including [8; 9; 10] for image generation. In this paper, we introduce interferometric NNs (INNs) consisting of interferometers in Sec. II and then interferometric GANs (IGANs) composed of INNs in Sec. III.

Interferometers play a pivotal role in metrology, astronomy, optics, and quantum physics. Feynman, in his renowned lectures, introduced quantum (wave-particle) behavior through double-slit (two-path) interference experiments [11]. The wave-particle duality, a subject of the famous Bohr-Einstein debates, has been further investigated by numerous scientists using two- and multi-path interferometers (for instance, see [12; 13; 14; 15; 16; 17; 18; 19]). Our INNs are made of sequences of such multi-path interferometers, where each interferometer is a combination of beamsplitters (or fiber couplers) and a phase shifter [14; 17; 18; 20]. A beamsplitter is represented by a discrete Fourier transformation, and the phase shifters hold all the learnable parameters of an INN. In principle, an INN can work with classical or quantum waves. However, a sequence of interferometers can be viewed as a parameterized quantum circuit (PQC), where a Fourier transformation can be implemented exponentially faster than its classical counterpart thanks to Peter Shor [21].

The PQCs [22] are integral components of the variational quantum algorithms (VQAs), including the variational quantum eigensolver (VQE) [23; 24] and its variants [25; 26; 27; 28] for determining molecular energies and the quantum approximate optimization algorithm (QAOA) [29] for solving combinatorial optimization problems by recasting them as energy minimization in Ising spin glass systems [30]. Recognizing that fully fault-tolerant quantum computers may still be decades away, VQAs are harnessing the capabilities of today's noisy intermediate-scale quantum (NISQ) computers [31] by providing a framework to build a wide range of applications. Some of these applications include finding ground [23; 24] and excited [32; 33] energy states, combinatorial optimization [29; 30], simulating quantum dynamics [34; 35; 36; 37; 38], solving a system of linear [39; 40; 41] and non-linear [42; 43] equations, data classification [44; 45; 46; 47; 48; 49; 50; 51], and generative tasks [52; 53; 54; 55; 56; 57]. For an in-depth review of VQAs and NISQ algorithms, we refer readers to [68; 69]. An INN also belongs to the VQA family, as it can provide an ansatz for a VQA.

We employ INNs for combinatorial optimization and image classification in Sec. II and IGANs for image generation in Sec. III. The notion of quantum GANs was initially introduced in [53] and numerically implemented in [54] using PQCs. Subsequently, various quantum GANs have been proposed to generate samples from random distributions [55; 56; 57], as well as to generate images [58; 59; 60; 61; 62]. Leveraging the Born rule, which provides a probability distribution from a (parameterized) quantum state, PQCs as Born machines [63; 64; 65; 66] are utilized for these generative tasks.

We provide a summary of our contributions at the beginning of the next two sections and an overarching conclusion and outlook in Sec. IV. One can find a dedicated Jupyter Notebook for each of the problems discussed in this paper within our GitHub repository [70]. These problems are tackled using INNs, which are classically simulated in the PyTorch machine-learning framework [71]. PyTorch automatically manages the gradient of a loss or an energy function. For the training of each INN, we employ the Adam optimizer [72].
Inn

In this section, our primary contributions are in (3)-(6), (14), Figs. 1-5, Table 1, and [70]. Here we introduce INNs through (3)-(6), (14), Figs. 1, 4, and 5. Initially applied to specific combinatorial optimization problems, we illustrate their performance in Fig. 2. Subsequently, we use them for image classification tasks, and their performance is detailed in Fig. 3 and Table 1.

The INN in Fig. 1 is made of a sequence of interferometers, where each interferometer has \(d\) distinct paths represented by the computation (orthonormal) basis

\[\mathcal{B}=\left\{\left|0\right\rangle,\cdots,\left|d-1\right\rangle\right\} \tag{1}\]

of the associated Hilbert space. Basically, this whole setup is a \(d\)-level quantum computer. In the case of \(d=2^{n}\), one can construct an equivalent PQC of \(n\) qubits (2-level quantum systems) as shown in Fig. 1.

In a machine-learning task, one has to feed the data to a NN, and the data is a collection of feature vectors \(\mathbf{x}=(x_{0},\cdots,x_{d-1})\in\mathbb{R}^{d}\). Whereas, in the case of INN, we have to encode every \(\mathbf{x}\) into an input state vector \(\left|\mathbf{x}\right\rangle\) through a unitary transformation. Such a process is called feature encoding [51]. Throughout the paper, we consider the _amplitude_ encoding

\[\left|\mathbf{x}\right\rangle=\frac{1}{\left|\mathbf{x}\right|}\sum_{j=0}^{d- 1}x_{j}\left|j\right\rangle, \tag{2}\]

where \(\left|\mathbf{x}\right|\) is the Euclidean norm of \(\mathbf{x}\).

As shown in Fig. 1 following [14; 17; 18; 20], an interferometer has two essential elements: _beamsplitters_ (or fiber couplers) and a _phase shifter_, which, in our case, are mathematically described by the unitary operators

\[F =\frac{1}{\sqrt{d}}\sum_{k=0}^{d-1}\sum_{j=0}^{d-1}\omega^{k\,kj} \left|k\right\rangle\!\!\left\langle j\right|\quad\text{and}\] \[U(\boldsymbol{\phi}) =\sum_{j=0}^{d-1}\mathrm{e}^{\mathrm{i}\phi_{j}}|j\rangle\! \!\left\langle j\right|, \tag{3}\]

respectively. In (3), \(\omega=\mathrm{e}^{\mathrm{i}\frac{2\pi}{d}}\) is a \(d\)th root of unity, \(\mathtt{i}=\sqrt{-1}\), and \(\boldsymbol{\phi}=(\phi_{0},\cdots,\phi_{d-1})\in\mathbb{R}^{d}\) carry phases. Here the discrete Fourier transformation \(F\) acts as a generalization of a symmetric 50/50 beamsplitter for a 2-path interferometer. It turns a single path \(|j\rangle\) into an equal superposition of all the paths of \(\mathcal{B}\). The phase shifter \(U(\boldsymbol{\phi})\) shifts the phase by \(\phi_{j}\) in the \(j\)th path. By the way, using the operators from (3), one can realize a phase encoding \(\left|\mathbf{x}\right\rangle_{\mathrm{ph}}=U(\mathbf{x})F|0\rangle\) of the data, which might be useful in some applications. We want to emphasize that two layers of the Hadamard and diagonal unitary operations like \(U\) are used in [48] to present a feature map that is hard to simulate classically. Additionally, parameterized diagonal unitary operations such as \(U(\boldsymbol{\phi})\) have been employed to make a quantum simulation faster in [37; 38] (see also [36]).

For every problem, we have executed \(F\) in [70], utilizing the fast Fourier transform available in PyTorch. However, its implementation can reach exponential speedup on a quantum computer due to Shor's algorithm for quantum Fourier transforms [21]. The phase shifter can be decomposed as a product of controlled-phase gates [36]:

\[U(\boldsymbol{\phi}) =\prod_{j=0}^{d-1}\mathrm{cp}(\phi_{j})\,,\quad\text{where}\] \[\mathrm{cp}(\phi_{j}) =\mathrm{e}^{\mathrm{i}\phi_{j}}|j\rangle\!\!\left\langle j \right|+\left(I-|j\rangle\!\!\left\langle j\right|\right)=\mathrm{e}^{ \mathrm{i}\phi_{j}|j\rangle\!\!\left\langle j\right|} \tag{4}\]

and \(I\) is the identity operator. In the case of \(d=2^{n}\), with the binary representation of integers \(j=\sum_{\mu=1}^{n}2^{n-\mu}j_{\mu}\) and association \(|j\rangle\equiv\otimes_{\mu=1}^{n}|j_{\mu}\rangle\), one can show that \(\mathrm{cp}\) is indeed a controlled-phase gate with \(n-1\) control qubits. Although there are exponentially many \(\mathrm{cp}\) gates relative to \(n\) in the decomposition, they can all, in principle, be simultaneously realized due to their commutativity. By the way, for \(\phi_{j}=\pi\), \(\mathrm{cp}(\phi_{j})\) becomes a quantum oracle of Grover's algorithm [73].

Now, we have all the elements to describe the INN of Fig. 1. The INN is made of \(L\) layers of interferometers, where the output from one is fed to the next. An input to the first interferometer is \(\left|\mathbf{x}\right\rangle\) from (2), the output from

Figure 1: The top and bottom pictures show an INN and its equivalent PQC in the case of \(d=2^{n}\), respectively. The INN is a sequence of \(L\) interferometers, where the output from one is fed to the next. Each interferometer has \(d\) distinct paths colored in green, beamsplitters (or fiber couplers) exhibited by blue-colored slabs, and a phase shifter shown by red disks in the top picture. In the bottom picture, brown wires depict \(n\) qubits, and the beamsplitter and phase shifter are displayed by blue and red-rounded rectangles, respectively. An input to the INN or to the PQC is a ket \(\left|\mathbf{x}\right\rangle\) from (2), the beamsplitter \(F\) and phase shifters \(U(\boldsymbol{\phi}^{\downarrow})\) are given in (3), and the output \(\left|\mathbf{x},\Phi\right\rangle\) is read by the detectors shown in black. The whole setup is summarized in (5).
the last one is

\[\ket{\mathbf{x},\Phi} =\left(\prod_{l=1}^{L}FU(\mathbf{\phi}^{l})\right)F\ket{\mathbf{x}}, \quad\text{and}\] \[p =\ket{\langle 0|\mathbf{x},\Phi\rangle|^{2}} \tag{5}\]

is the probability of getting 0th outcome (that is, detecting quantum particles in the 0th path). Equivalently, in the case of \(n\) qubits, we have to perform measurements on all the qubits to get \(p\) as \(\ket{0}=\ket{0\cdots 0}\) in the binary representation. Equations in (5) completely specify the INN in Fig. 1, and all its learnable parameters are denoted by \(\Phi=(\mathbf{\phi}^{1},\cdots,\mathbf{\phi}^{L})\in\mathbb{R}^{L\times d}\), which carries \(L\) phase-vectors, one for each interferometer.

One can view the INN as a differentiable function \(p:\mathbb{R}^{L\times d}\rightarrow[0,1]\) on the parameter space. In optimization or machine-learning tasks, PyTorch automatically manages the gradient computation of an energy or a loss function [given in (8), (11), and (17)] of \(p\) with respect to parameter updates. Without such a framework, we would need to compute derivatives such as

\[\frac{\partial p}{\partial\phi_{j}^{l}} =\bra{\mathbf{x}}\cdots\partial_{j}U^{\dagger}(\mathbf{\phi}^{l}) \cdots\ket{0}\!\bra{0}\cdots U(\mathbf{\phi}^{l})\cdots\ket{\mathbf{x}}+\] \[\bra{\mathbf{x}}\cdots U^{\dagger}(\mathbf{\phi}^{l})\cdots\ket{0} \!\bra{0}\cdots\partial_{j}U(\mathbf{\phi}^{l})\cdots\ket{\mathbf{x}},\text{ where}\] \[\partial_{j}U(\mathbf{\phi}^{l})=\frac{\partial U}{\partial\phi_{j}^{ l}}=\mathbbm{i}\,\mathrm{e}^{i\phi_{j}^{l}}|j\rangle\!\langle j|=\frac{\partial \mathrm{CP}(\phi_{j}^{l})}{\partial\phi_{j}^{l}}\,. \tag{6}\]

The second equation in (6) expresses the derivative of the \(l\)th phase shifter with respect to the phase \(\phi_{j}^{l}\), and it also turns out to be the derivative of the cp operator. Similarly, one can get the derivative \(\partial_{j}U^{\dagger}(\mathbf{\phi}^{l})\) of the adjoint (denoted by \(\dagger\)) of \(U\).

Now, we demonstrate applications of the INN first for quadratic unconstrained binary optimization (QUBO) and then for classification. A \(n\)-variable QUBO problem is equivalent to finding a minimum energy eigenstate of an associated \(n\)-qubit Hamiltonian (for many such problems, see [30])

\[H=\sum_{\mu,\nu=1}^{n}\mathsf{Q}_{\mu\nu}P_{\mu}P_{\nu}\,, \tag{7}\]

where \(\mathsf{Q}\) is a \(n\times n\) real matrix that defines the problem, and the projector \(P_{\mu}\) acts on the \(\mu\)th qubit as \(P_{\mu}|j_{\mu}\rangle=j_{\mu}|j_{\mu}\rangle,j_{\mu}\in\{0,1\}\). By the binary representation \(j=\sum_{\mu}2^{n-\mu}j_{\mu}\) and \(|j\rangle=\otimes_{\mu}|j_{\mu}\rangle\), \(H\) is diagonal in the computational basis \(\mathcal{B}\) of (1). The global minimum will be \(E_{j}=\langle j|H|j\rangle\) for some \(j\in\{0,\cdots,d-1\}\). Since \(\{E_{j}\}\) does not have any structure like convex functions, finding the minimum out of exponentially (\(d=2^{n}\)) many possibilities is an NP-hard problem.

Several successful VQAs, including the VQE [23; 24] and QAOA [29], have been proposed to find good approximate solutions. Their basic structure involves preparing a parametric state (referred as ansatz), computing and minimizing the energy expectation value, and obtaining an approximate solution from the optimized ansatz. A QAOA ansatz is constructed by sequentially applying the problem and a mixer unitary operators, whereas the VQE ansatz is generated by applying unitary coupled cluster operators to an initialized state (for more details, see [68; 24; 69]). Our approach is similar as narrated next, with the key distinction being that our ansatz is derived from the INN of Fig. 1.

We start with the input ket \(\ket{\mathbf{x}}=\ket{0}\) and reach the output ket \(\ket{0,\Phi}\equiv\ket{\Phi}\) as per (5). Then we perform measurements to compute the energy expectation value \(\mathcal{E}(\Phi)\) to minimize it over the parameter space and obtain

\[\Phi_{\text{sol}} =\underset{\Phi}{\text{argmin}}\ \mathcal{E}(\Phi)\,, \tag{8}\] \[\mathcal{E}(\Phi) =\bra{\Phi}H|\Phi\rangle=\sum_{j=0}^{d-1}E_{j}|\langle j|\Phi \rangle|^{2}=\sum_{\mu,\nu=1}^{n}\mathsf{Q}_{\mu\nu}\langle P_{\mu}P_{\nu} \rangle\,.\]

Usually, quantum and classical computers are put in a loop to compute and minimize \(\mathcal{E}\), respectively. To compute \(\mathcal{E}\), one needs only \(n^{2}\) expectation values \(\langle P_{\mu}P_{\nu}\rangle\) of local observables. After (8), we perform measurements on the optimized ansatz \(\ket{\Phi_{\text{sol}}}\) in the computational basis \(\mathcal{B}\) to get the most probable outcome

\[j_{\text{sol}}=\underset{j}{\text{argmax}}\ \lvert\langle j|\Phi_{\text{sol}} \rangle\rvert^{2} \tag{9}\]

as our solution. If the global minimum \(E_{\text{min}}=\min_{j}\{E_{j}\}\) is known then one can compute

\[\text{optimality gap}=\left|\frac{E_{j_{\text{sol}}}-E_{\text{min}}}{E_{ \text{min}}}\right|\times 100\,, \tag{10}\]

where \(E_{j_{\text{sol}}}=\langle j_{\text{sol}}|H|j_{\text{sol}}\rangle\) corresponds to the obtained solution via the INN.

For \(n=17\) qubits, to generate QUBO instances, we sampled \(\mathsf{Q}\) from one of the two distributions: the uniform

Figure 2: The plot illustrates optimality gaps for two distinct sets of QUBO problem instances, each containing one thousand instances for \(n=17\) qubits. The red circles (\(\circ\)) and blue triangles (\(\blacktriangle\)) represent the gaps for the two sets where the \(\mathsf{Q}\) matrix is sampled from the discrete uniform and standard normal distributions, respectively. The INN consistently yields the global minimum, achieving a success rate of 95% and 89% for the uniform and normal distributions, respectively. In both cases, the optimality gaps remain under 2%.
distribution over the set \(\{-10,\cdots,10\}\) of integers and the standard normal distribution. From each distribution, we generated one thousand QUBO instances. For each instance, we obtained \(E_{\mathrm{min}}\) through an exhaustive search over \(d=2^{17}=131072\) possibilities and reported the optimality gap in Fig. 2. In the case of discrete uniform distribution, we achieved the global minimum in 958 out of 1000 instances, meaning that the gap was \(\leqslant 10^{-10}\). Over all the instances, the mean and maximum gaps were 0.04% and 1.72%, respectively. When using the standard normal distribution, we reached the global minimum in 894 instances, and the mean and maximum gaps were 0.04% and 1.52%, respectively. For each instance, we maintained a consistent configuration with a number of layers, \(L=2\), and carried out 201 epochs of training with a learning rate of 0.05 and \(\mathrm{betas}=(0.5,0.9)\) for the Adam optimizer. Further details can be found in [70].

It will be interesting to see how the INN will perform if we increase the system size \(n\). For this, we need _quantum_ hardware that can store \(|\Phi\rangle\) and return \(\mathcal{E}(\Phi)\). Beyond certain \(n\), classical hardware cannot store exponentially many numbers \(\langle j|\Phi\rangle\) for \(j=0,\cdots,2^{n}-1\). It is a true power for a quantum computer [23]. We are not exploring this direction any further but moving to the next problem.

In a \(C\)-class classification problem, our goal is to predict the true label \(y=(y_{1},\cdots,y_{C})\in\{0,1\}^{C}\) for a data point based on its features provided in \(\mathbf{x}\). When the true class is \(c\), \(y_{c}=1\), and the remaining components of \(y\) are set to zero. During the training process of a NN for the problem, the objective is to minimize a specific loss function, such as the cross-entropy (negative log-likelihood)

\[\mathcal{L}=-\frac{1}{N}\sum_{i=1}^{N}\sum_{c=1}^{C}y_{c}^{[i]}\ln(p_{c}^{[i]}) \tag{11}\]

within the parameter space. The outer summation in (11) is performed over a mini-batch \(\{\mathbf{x}^{[i]},y^{[i]}\}_{i=1}^{N}\) containing \(N\) data points.

Suppose the dataset contains only \(C=2\) classes, then one can interpret \(p\) of (5) as the probability of a data point belonging to the positive class (\(c=1\)) and \(1-p\) as the probability of it belonging to the negative class (\(c=0\)). Consequently, the INN exhibited in Fig. 1 can be employed for binary classification. To illustrate this, we have utilized the MNIST dataset [74], which contains images of 0 to 9 handwritten digits. To create a 2-class classification problem, we gathered all the images of only two specific digits, denoted as \(a\) and \(b\), which are respectively labeled as the negative class \(y=(1,0)\) and the positive class \(y=(0,1)\). After the training of the INN, we evaluate its performance by using

\[\mathrm{accuracy} = \frac{\mathrm{TP}+\mathrm{TN}}{\mathrm{TP}+\mathrm{TN}+\mathrm{ FP}+\mathrm{FN}}\quad\mathrm{and}\] \[\mathrm{f}_{1} = \frac{2\,\mathrm{TP}}{2\,\mathrm{TP}+\mathrm{FP}+\mathrm{FN}} \tag{12}\]

on the test set. The results are presented in Fig. 3 for every \(a,b\in\{0,\cdots,9\}\) provided \(a\neq b\). Here, TP, TN, FP, and FN represent the counts of true positives, true negatives, false positives, and false negatives, respectively. Note that if we were to interchange the labels of \(a\) and \(b\), then \(\mathcal{L}\), the trained INN, accuracy, and \(\mathrm{f}_{1}\) score will be different. As a result, the matrices in Fig. 3 are not symmetric around their diagonals. Various binary classifications have been carried out on the MNIST dataset in [44; 45; 46; 49; 51] employing diverse PQCs. In our case, we achieve comparable performance, as illustrated in Fig. 3.

An image consists of \(d=\mathsf{C}\times\mathsf{H}\times\mathsf{W}\) pixel values, where \(\mathsf{C}\), \(\mathsf{H}\), and \(\mathsf{W}\) denote the number of channels, height, and width in pixels, respectively. Since the MNIST dataset contains black and white images, \(\mathsf{C}=1\). To get the results of Fig. 3, we have pre-processed the data, setting \(\mathsf{H}=\mathsf{W}=2^{4}\), and scaled all pixel values to the range of \(-1\) to \(1\). After flattening an image, we obtain a \(d=2^{8}\) component feature vector \(\mathbf{x}\), which is then transformed into

Figure 3: The top and bottom matrices display the accuracy and \(\mathrm{f}_{1}\) scores of trained INNs for binary classification problems. The entry in \(a\)th row and \(b\)th column corresponds to the image classification task of digit \(a\) (negative class) against digit \(b\) (positive class). It can be observed that the INN achieves accuracy ranging from 0.81 to 0.99 and \(\mathrm{f}_{1}\) scores between 0.75 and 0.99.
\(|\mathbf{x}\rangle\) following (2). Subsequently, we pass \(|\mathbf{x}\rangle\) through the INN shown in Fig. 1, utilizing \(L=2\) layers. Other hyperparameters, including the learning rate, betas, batch size, and number of epochs, are set to \(0.01\), \((0.5,0.9)\), \(2^{6}\), and \(10\), respectively (for additional details, refer to [70]). For every \((a,b)\), we maintain the same settings as described above to obtain the results presented in Fig. 3. In every case, the INN achieves an accuracy of more than \(80\%\) and \(\mathrm{f}_{1}\) score more than \(75\%\). These scores can potentially be further improved by modifying the INN architecture and fine-tuning the hyperparameters.

To classify images of the MNIST dataset into the \(C=10\) classes, we need a NN that takes a \(d_{\mathrm{in}}\)-dimensional input and provides a \(d_{\mathrm{out}}\)-dimensional output, where \(d_{\mathrm{out}}\) is not necessarily the same as \(d_{\mathrm{in}}\). Let us take \(d_{\mathrm{in}}=3\), \(d_{\mathrm{out}}=2\), and compare

\[\mathbf{x}^{\prime}=W\mathbf{x}=\begin{pmatrix}w_{00}&w_{01}&w_{02}\\ w_{10}&w_{11}&w_{12}\end{pmatrix}\begin{pmatrix}x_{0}\\ x_{1}\\ x_{2}\end{pmatrix}\quad\text{with} \tag{13}\]

The above equations represent a single linear layer of a classical NN (multi-layer perceptron) [1, 2, 3] and a single layer of INN of Fig. 1, respectively. Here, we have used amplitude encoding (2) by assuming \(\|\mathbf{x}\|=1\). The linear layer combines the input features from \(\mathbf{x}\) in a weighted manner to generate new features in \(\mathbf{x}^{\prime}\), while the INN layer creates _interference_ among the input features to generate new features.

Two more differences can be observed through (13). Firstly, the output dimension is \(2\) in the case of the linear layer, while it is \(3\) (not equal to the desired \(d_{\mathrm{out}}\)) for the INN layer. Secondly, each row of \(W\) is _unrelated_ to the others, whereas the rows of any unitary matrix must follow the _orthogonality relation_. As a result, we get _distinct_ new features \(\sum_{j}w_{0j}x_{j}\) and \(\sum_{j}w_{1j}x_{j}\) from the linear layer. Whereas, from the INN layer, the new features like \(\sum_{j}\mathrm{e}^{\mathrm{i}\phi_{j}}x_{j}\) and \(\sum_{j}\mathrm{e}^{\mathrm{i}\phi_{j}}x_{j}\,\omega^{j}\) represent very similar functions (machine learning modes) of the _same_ set \(\mathbf{\phi}\) of parameters, and \(\omega=\mathrm{e}^{\mathrm{i}\frac{2\pi}{3}}\) is a constant. Furthermore, their derivatives with respect to a parameter only differ by a constant. So, to get \(d_{\mathrm{out}}\)_distinct_ features, we put \(d_{\mathrm{out}}\) distinct INNs of Fig. 1 in _parallel_ and obtain

\[|\mathbf{x},\Phi^{t}\rangle =\mathbf{U}(\Phi^{t})|\mathbf{x}\rangle\,,\quad\text{where}\] \[\mathbf{U}(\Phi^{t}) =\left(\prod_{l=1}^{L}FU(\mathbf{\phi}^{t_{l}})\right)F\quad\text{and}\] \[p_{t} =|\langle 0|\mathbf{x},\Phi^{t}\rangle|^{2} \tag{14}\]

for \(t=0,\cdots,d_{\mathrm{out}}-1\). Similar to how (5) represents a sequence of interferometers, (14) portrays a block of sequences of interferometers, as displayed in Fig. 4. In the figure, we also illustrate the similarity of (14) to a linear layer that represents \(d_{\mathrm{out}}\) linear regression models in _parallel_. One can observe that, from the same \(|\mathbf{x}\rangle\), we get \(d_{\mathrm{out}}\) distinct probabilities \(p_{t}\) through \(d_{\mathrm{out}}\)_mutually independent_ unitary operators \(\mathbf{U}(\Phi^{t})\), each of which represents a sequence of \(L\) interferometers. Since these probabilities are independent of each other, they are not required to sum up to \(1\). If their normalization is needed, one can apply the function

\[\mathrm{softmax}\left(p_{t}\right)=\frac{\mathrm{e}\,p_{t}}{\sum_{t=0}^{d_{ \mathrm{out}}-1}\mathrm{e}\,p_{t}} \tag{15}\]

to every component of \(\mathbf{p}=(p_{0},\cdots,p_{d_{\mathrm{out}}-1})\). Essentially, the INN block (14) creates distinct new features \(\mathbf{p}\) from the old \(\mathbf{x}\) through independent interferences. The vector \(\mathbf{p}\) can serve as a final output or an input to the next INN block as shown in Fig. 5.

Figure 5 shows an INN whose architecture is defined by \(\mathbf{d}=(d_{0},\cdots,d_{M-1},d_{M})\) and \(\mathbf{L}=(L_{0},\cdots,L_{M-1})\), where the number of blocks \(M\) is \(2\). Its \(m\)th block takes \(d_{m}\)-dimensional input and gives \(d_{m+1}\)-dimensional output. Within the block, there are \(d_{m+1}\) parallel sequences of interferometers represented by \(\otimes_{t=0}^{d_{m+1}-1}\mathbf{U}(\Phi^{t})\), and the length of each sequence is \(L_{m}\).

While the INN in Fig. 5 draws inspiration from multilayer perceptrons (NNs) [1, 2, 3] like the quantum NNs (QNNs) in [75, 76, 44], it differs in the following aspects. A single quantum preceptron in [75] is a sum of single-qubit unitary operators, whereas it is a product of unitary operators in the case of INN as shown in (5). Compared

Figure 4: The top and bottom pictures exhibit a linear layer [for example, given in (13)] of a classical NN and an interferometric block [specified by (14)] of an INN, respectively. By the color coding one can see a similar _parallel_ structure in both cases. The linear layer and INN block turn the input \(\mathbf{x}\in\mathbb{R}^{d_{\mathrm{in}}}\) into the outputs \(\mathbf{x}^{\prime}\in\mathbb{R}^{d_{\mathrm{out}}}\) and \(\mathbf{p}\in[0,1]^{\times d_{\mathrm{out}}}\), respectively. In the bottom picture, the squares represent \(d_{\mathrm{in}}\times d_{\mathrm{in}}\)_independent_ unitary matrices \(\mathbf{U}(\Phi^{t})\), each of which is a sequence of interferometers as depicted in Fig. 1.
to [44; 76], the INN has a parallel structure of unitary operators, as depicted in Fig. 4. In contrast to the QNNs from [75; 76; 44], where each node represents a qubit, INNs do not necessitate a qubit system. INNs are applicable in any dimension and explicitly incorporate Fourier transformations.

We have employed INNs, as illustrated in Fig. 5, for image classification on both the MNIST and FashionMNIST [77] datasets, each with \(C=10\) classes. The FashionMNIST dataset, like MNIST, consists of images representing clothing items from ten (labeled \(0\) to \(9\)) different categories. In both cases, we have adopted \(\mathbf{d}=(2^{8},2^{6},C)\), the amplitude encoding of \(16\times 16\) images, \(0.01\) learning rate, \(2^{5}\) batch size, and \(3\) epochs (for further details, refer to [70]). Their performances on the test sets are presented in Table 1 for \(\mathbf{L}=(1,1)\) and \((2,2)\). As we increase the number of layers, the performance--measured by the accuracy and average f\({}_{1}\)--of the INN improves for both datasets. In summary, we have attained accuracies and average f\({}_{1}\) scores of \(93\%\) and \(83\%\) on the MNIST and FashionMNIST datasets, respectively.

The f\({}_{1}\) score in (12) pertains to the positive class. We can utilize a similar formula to calculate the f\({}_{1}\) score for each class and then derive the average f\({}_{1}\) score. It is also possible to extend the accuracy formula presented in (12) from \(C=2\) to \(C=10\) classes. In the case of the cross-entropy loss function described in (11), it is essential to ensure that the probabilities are appropriately normalized. This can be accomplished by employing the softmax function outlined in (15).

We conclude this section with the following two observations: (i) It is of interest to investigate whether the pair \(\{U(\mathbf{\phi}),F\}\) from (3) constitutes a universal set--capable of generating any unitary operation through multiplications--for quantum computation with a \(d\)-level system. Notably, for \(d=2\), it is a known universal set [78]. Moreover, \(U(\mathbf{\phi})\) encompasses all diagonal unitary operations and can be expressed as a linear combination of different powers of \(Z=\sum_{j=0}^{d-1}\omega^{j}|j\rangle\langle j|\). Through multiplication, \(Z\) and \(X=F^{3}ZF\) can generate the Heisenberg-Weyl group [79], whose elements constitute the unitary operator bases [80].

(ii) As a black and white image is a two-dimensional (2D) object, instead of (2), one can use

\[|\mathbf{x}\rangle_{\text{2D}}=\frac{1}{\|\mathbf{x}\|}\sum_{j=0}^{\mathsf{H- 1}}\sum_{j^{\prime}=0}^{\mathsf{W-1}}x_{jj^{\prime}}\,|j\rangle\otimes|j^{ \prime}\rangle\,, \tag{16}\]

to encode an image \(\mathbf{x}\) into two quantum subsystems with dimensions \(\mathsf{H}\) and \(\mathsf{W}\), where \(d=1\times\mathsf{H}\times\mathsf{W}\). Then, one can use the tensor product \(F_{\mathsf{H}}\otimes F_{\mathsf{W}}\) and \(U(\mathbf{\phi})=\sum_{k}\sum_{k^{\prime}}\mathrm{e}^{\mathrm{i}\phi_{kk^{\prime }}}|k\rangle\!\langle k|\otimes|k^{\prime}\rangle\!\langle k^{\prime}|\) at the place of (3). The tensor product performs a 2D discrete Fourier transform on the matrix \(x_{jj^{\prime}}\), with \(F_{\mathsf{H}}\) and \(F_{\mathsf{W}}\) operating exclusively on their respective subsystems. Afterward, \(U(\mathbf{\phi})\) executes a componentwise multiplication between the phase matrix (filter) \(\mathrm{e}^{\mathrm{i}\phi_{kk^{\prime}}}\) and the Fourier-transformed matrix \(\widehat{x}_{kk^{\prime}}\). As per the convolution theorem, such multiplication in the frequency domain is linked to _convolution_ in the spatial domain through Fourier transforms. Fourier transforms are fundamental components of digital image processing, particularly for designing various filters in the frequency domain [81].

## III Igan

In this section, we introduce IGANs, and our primary contributions are summarized as follows: Tables 2 and 3 detail the architectures and hyperparameters of IGANs, Algorithm 1 outlines the training process, and the implementation can be found in [70]. Figure 9 offers insights into how losses and probabilities evolve during the training. And, Figs. 7 and 8 display sample images generated after the training process.
generated data and provides a probability score for an input being real.

To enhance classification accuracy, the discriminator aims to drive \(D(\mathbf{x})\) closer to 1 through maximization and \(D(G(\mathbf{z}))\) closer to 0 through minimization. Both objectives are accomplished by minimizing the discriminator's loss

\[\mathcal{L}_{D} =-\frac{1}{2N}\sum_{i=1}^{N}\left[\ln(D(\mathbf{x}^{[i]}))+\ln(1-D (G(\mathbf{z}^{[i]})))\right]\,,\] \[\mathcal{L}_{G} =-\frac{1}{N}\sum_{i=1}^{N}\ln(D(G(\mathbf{z}^{[i]}))) \tag{17}\]

serves as a guide for the generator. It encourages the generator to produce increasingly realistic data, thereby pushing \(D(G(\mathbf{z}))\) closer to 1. In (17), every summation is over a mini-batch, where \(\{\mathbf{x}^{[i]}\}_{i=1}^{N}\) originates from the given data and \(\{\mathbf{z}^{[i]}\}_{i=1}^{N}\) is drawn from a prior (in our case, the standard normal) distribution.

The two networks are trained simultaneously as described in Algorithm 1 inspired from [7]. Initially, \(\mathcal{L}_{D}\) is minimized with respect to the discriminator's parameters, followed by minimizing \(\mathcal{L}_{G}\) with respect to the generator's parameters within a single iteration (epoch). The two networks compete with each other and seek the Nash equilibrium point, where the discriminator is unable to distinguish between real and fake data. There both the probabilities \(D(\mathbf{x})\) and \(D(G(\mathbf{z}))\) attain a value of \(\frac{1}{2}\), resulting in \(\mathcal{L}_{D}=\ln(2)=\mathcal{L}_{G}\).

Training a GAN can be a challenging and unstable process. Therefore, several GAN variants, such as Deep Convolutional GANs [8] and Wasserstein GANs [9; 10], have been developed. Deep Convolutional GANs introduced architectural guidelines for generator and discriminator networks, leading to more stable training and the generation of realistic and detailed images. Wasserstein GANs introduced the Wasserstein distance as a more stable and meaningful loss function for GAN training, effectively addressing issues like mode collapse.

```
Input: batch size \(N\), number of epochs, learning rate, \(\mathrm{Adam}\) hyperparameters \((\beta_{1},\beta_{2})\) from Table III. \(\bullet\) Initialize parameters of discriminator \(\boldsymbol{\Phi}_{D}\) and of generator \(\boldsymbol{\Phi}_{G}\) with the Xavier initialization [82]. fornumber of epochsdo fornumber of mini-batchesdo\(\bullet\) Sample mini-batch \(\{\mathbf{z}^{[i]}\}_{i=1}^{N}\) from the standard normal distribution \(\mathcal{N}(0,1)\). \(\bullet\) Sample mini-batch \(\{\mathbf{x}^{[i]}\}_{i=1}^{N}\) from the training dataset. \(\bullet\) Update \(\boldsymbol{\Phi}_{D}\leftarrow\mathrm{Adam}(\nabla_{\boldsymbol{\Phi}_{D}} \mathcal{L}_{D},\beta_{1},\beta_{2})\)\(\bullet\) Update \(\boldsymbol{\Phi}_{G}\leftarrow\mathrm{Adam}(\nabla_{\boldsymbol{\Phi}_{G}} \mathcal{L}_{G},\beta_{1},\beta_{2})\)  end for  end for
```

**Algorithm 1**IGAN training algorithm

In the quantum domain, several models have been employed for image generation, including the quantum GANs [58; 59; 60; 61; 62], Born machines [63; 64; 66], matrix product states [83], and quantum variational autoencoder [84]. In quantum GANs, described in [58; 59], images are created from multiple patches generated in parallel by sub-generators from a product state, which encodes the components of \(\mathbf{z}\) into rotation angles. The loss function in [59] is inspired by Wasserstein GANs [10]. In contrast, quantum GANs in [61; 62] leverage quantum state fidelity-based loss functions and incorporate principal component analysis for input image compression. In the case of hybrid quantum-classical GANs discussed in [60], a specific remapping method is employed to enhance the quality of generated images. The generative models in [63; 64; 83] are built using datasets of binary images. In [66; 84; 85], quantum models are employed for producing noise vectors \(\mathbf{z}\). Similar to our IGAN, the MNIST dataset is utilized in all the generative models introduced in [58; 59; 60; 61; 62; 66; 83; 84; 85].

Now, we introduce our IGANs, where both the generator \(G\) and discriminator \(D\) are INNs as shown in Fig. 5. It is important to note that our approach does not involve principal component analysis for input compression, binary images for training, or the inclusion of any classical layers within our INNs. For input loading, we always use amplitude encoding (2) for classical vectors such as \(\mathbf{x}\), \(\mathbf{z}\), and \(\mathbf{p}\). An output from an interferometric block is always a classical vector \(\mathbf{p}\) as illustrated in Fig. 4. From the generator, we get \(\mathbf{p}\in[0,1]^{\times d}\) with an image size of \(d=\mathsf{C}\times\mathsf{H}\times\mathsf{W}\), which, after the componentwise transformation

\[\mathbf{p}\rightarrow\tanh(2\,\mathbf{p}-1)=G(\mathbf{z})\in[-1,1]^{\times d }\,, \tag{18}\]

is fed to the discriminator during training. The discriminator's output, \(D(\mathbf{x})\) or \(D(G(\mathbf{z}))\), represents the probability, \(p\in[0,1]\), of the input being real, where the corresponding \(\mathbf{p}=(1-p\,,p)\). The training process for our

\begin{table}
\begin{tabular}{c|c|c|c} \hline \hline dataset & INN & \(\mathbf{d}\) & \(\mathbf{L}\) \\ \hline MNIST & \(D\) & \((20^{2},2^{6},1)\) & \((3,3)\) \\  & \(G\) & \((2^{3},2^{6},20^{2})\) & \((3,3)\) \\ \hline CelebA & \(D\) & \((32^{2},2^{7},1)\) & \((3,3)\) \\  & \(G\) & \((2^{5},2^{7},32^{2})\) & \((3,3)\) \\ \hline \hline \end{tabular}
\end{table}
Table II: The table presents the architecture of INNs for the discriminator \(D\) and generator \(G\) of separate IGANs trained on the two datasets. In both cases, every INN is of the structure depicted in Fig. 5.

\begin{table}
\begin{tabular}{c|c|c|c|c|c} \hline \hline dataset & data size & batch size & epochs & rate & \((\beta_{1},\beta_{2})\) \\ \hline MNIST & 5000 & \(2^{7}\) & 10 & 0.01 & \((0.5,0.9)\) \\ \hline CelebA & 50000 & \(2^{9}\) & 10 & 0.01 & \((0.5,0.9)\) \\ \hline \hline \end{tabular}
\end{table}
Table III: The table presents the hyperparameters utilized in Algorithm 1 for training separate IGANs on the two datasets.
IGANs is outlined in Algorithm 1 and its complete implementation is given in [70].

Separate IGANs were trained on the MNIST dataset [74] and CelebA datasets [86], which consists of images featuring celebrities' faces as displayed in Fig. 6. Subsequently, samples of generated images from their respective \(G\)s are presented in Fig. 7. In both cases, black and white images were used for training, and 2-block INNs were employed for both \(D\) and \(G\). For the two datasets, we have adopted slightly different architectures for the INNs, as detailed in Table 2, and slightly varied hyper-parameters, as specified in Table 3.

In the case of the MNIST dataset, we have randomly taken 5000 images of 0 to 9 digits for the training as specified in the "data size" column of Table 3. Each image has dimensions \(d=\mathsf{C}\times\mathsf{H}\times\mathsf{W}=1\times 20\times 20\), as evident from the first and last components of \(\mathsf{ds}\) for \(D\) and \(G\) in Table 2. In the top panel of generated images in Fig. 7, clear representations of numbers 0, 1, 3, 5, and 9 are easily noticeable. However, the images of 2, 4, and 6 appear less distinct, while those of 7 and 8 do

Figure 6: The top and bottom panels showcase samples of 64 _real_ images from the MNIST [74] and CelebA [86] datasets, respectively.

Figure 7: After separately training IGANs on the MNIST [74] and CelebA [86] datasets, the top and bottom panels showcase 64 _fake_ images generated by their respective generators. To enhance contrast, we have employed the transformation \(\mathbf{p}\rightarrow\tanh(3(2\,\mathbf{p}-1))\) at the place of (18) exclusively for image generation in the top panel and Fig. 8.
not seem to be present. In another set of experiments, we trained the same IGAN using images of a single digit at a time and achieved the results depicted in Fig. 8. These results demonstrate our capability to generate images corresponding to all ten digits, from 0 to 9. One can visually assess the quality of generated images by comparing them to real images through Figs. 6-8.

In the case of the CelebA dataset, we randomly selected 50000 images for training. These images were converted from color to black and white and resized to the dimensions of \(d=1\times 32\times 32\), as detailed in Table 2. With the table, one can also compute that the INNs for \(D\) and \(G\) have a total of 393600 and 405504 trainable parameters (phases), respectively. In Fig. 9, we depict how the discriminator and generator losses converge toward \(\ln(2)\) as training progresses. Additionally, we illustrate how the average discriminator probabilities

\[\mathcal{P}_{\text{real}}=\frac{1}{N}\sum_{i=1}^{N}D(\mathbf{x}^{[i]})\quad \text{and}\quad\mathcal{P}_{\text{fake}}=\frac{1}{N}\sum_{i=1}^{N}D(G(\mathbf{ z}^{[i]})) \tag{19}\]

approach \(\frac{1}{2}\). After the training, we got the ability to generate images featuring human faces, as showcased in Fig. (7). By the way, with the average probabilities, one can define \(\mathcal{L}^{\prime}_{D}=-\mathcal{P}_{\text{real}}+\mathcal{P}_{\text{fake}}\) and \(\mathcal{L}^{\prime}_{G}=-\mathcal{P}_{\text{fake}}\), inspired from the Wasserstein loss functions [9, 10, 59], at place of (17).

## IV Conclusion and outlook

We have introduced INNs, which are artificial neural networks composed of interferometers. An INN is a sequence of interferometric blocks, with each block containing parallel sequences of interferometers. While one could, in principle, use an INN with classical waves, it can be regarded as a quantum model, devoid of any classical layers.

We have shown that INNs are useful for optimization as well as supervised or unsupervised machine learning tasks. For the QUBO problems, we achieve the global

Figure 8: From top to bottom, each row displays eight _fake_ images generated by the \(G\) network trained separately on real images of digits \(0,\cdots,9\). For each digit, we maintained consistent architecture and hyperparameters (with the exception of data size) as those presented in Tables 2 and 3 for the MNIST dataset. The training data was sourced from the MNIST training set.

Figure 9: At the top, we present the evolution of \(\mathcal{L}_{D}\) and \(\mathcal{L}_{G}\) at various epochs, represented by orange and blue curves, respectively. Meanwhile, at the bottom, we illustrate the changes in \(\mathcal{P}_{\text{real}}\) and \(\mathcal{P}_{\text{fake}}\) over different epochs, denoted by red and green graphs, respectively. The definitions of these loss functions and average probabilities can be found in (17) and (19), respectively. In the top and bottom plots, the gray horizontal lines denote the values \(\ln(2)\) and \(\frac{1}{2}\), respectively. Both the plots pertain to the IGAN trained on the CelebA dataset.
minimum approximately 90% of the time, with the remaining instances typically falling within a range of 2% from the global minimum. In the context of multi-class image classification problems, we achieved accuracies of 93% and 83% on the MNIST and FashionMNIST datasets, respectively. While our accuracy falls short in comparison to state-of-the-art classical NNs, which provide 99.87% accuracy on the MNIST dataset [87], and 96.91% accuracy on the FashionMNIST dataset [88], it is important to note that INNs exhibit a simpler architecture. Nonetheless, our work represents a significant step forward in the development of more advanced quantum NNs.

We have introduced IGANs, made of INNs, for image generation. These IGANs have successfully generated images of 0 to 9 digits and human faces. While our image quality may currently lag behind that of classical GANs [7; 8; 9; 10], there is potential for enhancement through network modifications, architectural adjustments, and fine-tuning of hyperparameters. Last but not least, it is crucial to analyze the robustness of our INNs against noise and qubit errors, which will be the focus of future research.
