Classical and Quantum Algorithms for Orthogonal Neural Networks

The idea behind Orthogonal Neural Networks (OrthoNNs) is to add a constraint to the weight matrices corresponding to the layers of a neural network. Imposing orthogonality to these matrices has theoretical and practical benefits in the generalization error [3]. Orthogonality ensures a low weights redundancy and preserves the magnitude of the weight matrix’s eigenvalues to avoid vanishing gradients. In terms of complexity, for a single layer, the feedforward pass of an OrthoNN is simply a matrix multiplication, hence has a running time of $O(n^{2})$ if $n\times n$ is the size of the orthogonal matrix. It is also interesting to note that OrthoNNs have been generalized to Convolutional Neural Networks [4].

The main difficulty of OrthoNNs is to preserve the orthogonality of the matrices while updating them during gradient descent. Several algorithms have been proposed to this end [4, 6, 11], but they all point that pure orthogonality is computationally hard to conserve. Therefore, previous works allow for approximations: strict orthogonality is no longer required, and the matrices are often pushed toward orthogonality using regularization techniques during weights update.

We present two algorithms from [3] for updating orthogonal matrices.

The first algorithm is an approximate one, called Singular Value Bounding (SVB). It starts by applying the usual gradient descent update on the matrix, therefore making it not orthogonal anymore. Then, the singular values of the new matrix are extracted using Singular Value Decomposition (SVD), their values are manually pushed to be close to 1, and the matrix is recomposed hence enforcing orthogonality. This method shows less advantage on practical experiments [3]. It has a complexity of $O(n^{3})$ due to the SVD, and in practice is better than the second algorithm described below. Note that this running time is still asymptotically worse than $O(n^{2})$, the running time to perform standard gradient descent.

The second algorithm ensures perfect orthogonality by performing the gradient descent in the manifold of orthogonal matrices, called the Stiefel Manifold. In practice [3] this method showed a substantially advantageous classification results on standard datasets, indicating that perfect orthogonality can be a much more desired property than approximate one. This algorithm requires $O(n^{3})$ operations, and is quite prohibitive in practice. We give an informal step-by-step detail of this algorithm:

1. 1.

   Compute the gradient $G$ of the weight matrix $W$.

2. 2.

   Project the gradient matrix $G$ in the tangent space, (The space tangent to the manifold at this point $W$): multiply $G$ by some other matrices based on $W$:

   $$(I-WW^{T})G+\frac{1}{2}W(W^{T}G-G^{T}W)$$

   This requires several matrix-matrix multiplications. In the case of square $n\times n$ matrices, each has complexity $O(n^{3})$. the result of this projection is called the manifold gradient $\Omega$.

3. 3.

   update $W^{\prime}=W-\eta\Omega$, where $\eta$ is the chosen learning rate.

4. 4.

   Perform a retraction from the tangent space to the manifold. To do so we multiply $W^{\prime}$ by $Q$ factor of the QR decomposition, obtained using Gram Schmidt orthonormalization, which has complexity $O(2n^{3})$.

In conclusion, previous work points to the two following statements. First, perfect orthogonality can indeed increase classification accuracy. Second, previously known algorithms for training orthogonal neural networks were not as efficient as a standard training algorithm.

The quantum circuit proposed in this work (see Fig.1), which implements a fully connected neural network layer with an orthogonal weight matrix, uses only one type of quantum gate, the Reconfigurable Beam Splitter (RBS) gate. This two-qubit gate is parametrizable with one angle $\theta\in[0,2\pi]$. Its matrix representation is given as:

We can think of this gate as a rotation in the two-dimensional subspace spanned by the basis $\{\mathinner{|{01}\rangle},\mathinner{|{10}\rangle}\}$, while it acts as the identity in the remaining subspace $\{\mathinner{|{00}\rangle},\mathinner{|{11}\rangle}\}$. Or equivalently, starting with two qubits, one in the $\mathinner{|{0}\rangle}$ state and the other one in the state $\mathinner{|{1}\rangle}$, the qubits can be swapped or not in superposition. The qubit $\mathinner{|{1}\rangle}$ stays on its wire with amplitude $\cos\theta$ or switches with the other qubit with amplitude $+\sin\theta$ if the new wire is below ($\mathinner{|{10}\rangle}\mapsto\mathinner{|{01}\rangle}$) or $-\sin\theta$ if the new wire is above ($\mathinner{|{01}\rangle}\mapsto\mathinner{|{10}\rangle}$). Note that in the two other cases ($\mathinner{|{00}\rangle}$ and $\mathinner{|{11}\rangle}$) the $RBS$ gate acts as identity.

This gate can be implemented rather easily, either as a native gate, known as $FSIM$ [12], or using four Hadamard gates, two $R_{y}$ rotation gates, and two two-qubits CZ gates:

We now propose a quantum circuit that implements an orthogonal layer of a neural network. The circuit is a pyramidal structure of $RBS$ gates, each with an independent angle, as represented in Fig.3(a). In Section 2.3 and 3, more details are provided concerning respectively the input loading, and the equivalence with a neural network’s orthogonal layer.

To mimic a given classical layer with a quantum circuit, the number of output qubits should be the size of the classical layer’s output. We refer to the square case when the input and output sizes are equal, and to the rectangular case otherwise (Fig.4(a)).

The important property to note is that the number of parameters of the quantum pyramidal circuit corresponding to a neural network layer of size $n\times d$ is $(2n-1-d)*d/2$ exactly the same as the number of degrees of freedom of an orthogonal matrix of dimension $n\times d$.

For simplicity, we pursue our analysis using only the square case but everything can be easily extended to the rectangular case. As we said, the full pyramidal structure of the quantum circuit described above imposes the number of free parameters to be $N=n(n-1)/2$, the exact number of free parameters to specify a $n\times n$ orthogonal matrix.

In Section 3 we will show how the parameters of the gates of this pyramidal circuit can be easily related to the elements of the orthogonal matrix of size $n\times n$ that describes it. We note that alternative architectures can be imagined as long as the number of gate parameters is equal to the parameters of the orthogonal weight matrix and a simple mapping between them and the elements of the weight matrix can be found.

Note finally that this circuit has linear depth and is convenient for near term quantum hardware platforms with restricted connectivity. Indeed, the distribution of the $RBS$ gates requires only nearest neighbor connectivity between qubits.

Before applying the quantum pyramidal circuit, we will need to upload the classical data into the quantum circuit. We will use one qubit per feature of the data. For this, we use a unary amplitude encoding of the input data. Let’s consider an input sample $x=(x_{0},\cdots,x_{n-1})\in\mathbb{R}^{n}$, such that $\left\lVert x\right\rVert_{2}=1$. We will encode it in a superposition of unary states:

We can also rewrite the previous state as $\mathinner{|{x}\rangle}=\sum^{n-1}_{i=0}x_{i}\mathinner{|{e_{i}}\rangle}$, where $\mathinner{|{e_{i}}\rangle}$ represents the i^th unary state with a $\mathinner{|{1}\rangle}$ in the i^th position $\mathinner{|{0\cdots 010\cdots 0}\rangle}$. Recent work [13] proposed a logarithmic depth data loader circuit for loading such states. Here we will use a much simpler circuit. It is a linear depth cascade of $n$-1 $RBS$ gates which, due to the particular structure of our quantum pyramidal circuit, only adds 2 extra steps to our circuit.

The circuit starts in the all $\mathinner{|{0}\rangle}$ state and flips the first qubit using an $X$ gate, in order to obtain the unary state $\mathinner{|{10\cdots 0}\rangle}$ as shown in Fig.5. Then a cascade of $RBS$ gates allow to create the state $\mathinner{|{x}\rangle}$ using a set of $n-1$ angles $\alpha_{0},\cdots,\alpha_{n-2}$. Using Eq.(1), we will choose the angles such that, after the first $RBS$ gate of the loader, the qubits would be in the state $x_{0}\mathinner{|{100\cdots}\rangle}+\sin(\alpha_{0})\mathinner{|{010\cdots}\rangle}$ and after the second one in the state $x_{0}\mathinner{|{100\cdots}\rangle}+x_{1}\mathinner{|{010\cdots}\rangle}+\sin(\alpha_{0})\sin(\alpha_{1})\mathinner{|{001\cdots}\rangle}$ and so on, until obtaining $\mathinner{|{x}\rangle}$ as in Eq.(2). To this end, we simply perform a classical preprocessing to compute recursively the $n$-1 loading angles, in time $O(n)$. We choose $\alpha_{0}=\arccos(x_{0})$, $\alpha_{1}=\arccos(x_{1}\sin^{-1}(\alpha_{0}))$, $\alpha_{2}=\arccos(x_{2}\sin^{-1}(\alpha_{0})\sin^{-1}(\alpha_{1}))$ and so on.

The ability of loading data in such a way relies on the assumption that each input vector is normalized, i.e. $\left\lVert x\right\rVert_{2}=1$. This normalization constraint could seem arbitrary and impact the ability to learn from the data. In fact, in the case of an orthogonal neural network, this normalization shouldn’t degrade the training because orthogonal weight matrices are in fact orthonormal and thus norm-preserving. Hence, changing the norm of the input vector, by diving each component by $\left\lVert x\right\rVert_{2}$, in both classical and quantum settings is not a problem. The normalization would impose that each input has the same norm, or the same "luminosity" in the context of images, which can be helpful or harmful depending on the case.

It is important to notice that with our restriction to unary states, strong error mitigation techniques become available. Indeed, as we expect to obtain only quantum superposition of unary states at every layer, we can post process our measurements and discard the ones that present non unary states (i.e. states with more than one qubit in state $\mathinner{|{1}\rangle}$, or the ground state). The most expected error is a bit flip between $\mathinner{|{1}\rangle}$ and $\mathinner{|{0}\rangle}$. The case where two bit flips happen at the same time, which would change a unary state to a different unary state and would thus pass through our error mitigation, is even less probable. This error mitigation procedure can be applied efficiently to the results of a hardware demonstration and it has been used in the results presented in this paper.

As shown in Fig.8, when using the quantum circuit, the output is a quantum state $\mathinner{|{y}\rangle}=\mathinner{|{Wx}\rangle}$. As often in quantum machine learning [14], it is important to consider the cost of retrieving the classical outputs, using a procedure called tomography. In our case this is even more crucial since, between each layer, the quantum output will be converted into a classical one in order to apply a non linear function, and then reloaded for the next layer.

Retrieving the amplitudes of a quantum state comes at cost of multiple measurements, which requires running the circuit multiples times, hence adding a multiplicative overhead in the running time. A finite number of samples is also a source of approximation error in the final result. In this work, we will allow for $\ell_{\infty}$ errors [7]. The $\ell_{\infty}$ tomography on a quantum state $\mathinner{|{y}\rangle}$ with unary encoding on $n$ qubits requires $O(\log(n)/\delta^{2})$ measurements, where $\delta>0$ is the error threshold allowed. For each $j\in[n]$, $|y_{j}|$ will be obtained with an absolute error $\delta$, and if $|y_{j}|<\delta$, it will most probably not be measured, hence set to 0. In practice, one would perform as many measurements as is convenient during the experiment, and deduce the equivalent precision $\delta$ from the number of measurements made.

Note that it is important to obtain the amplitudes of the quantum state, which in our case are positive or negative real numbers, and not just the probabilities of the outcomes, which are the squares of the amplitudes. There are different ways of obtaining the sign of the amplitudes and we present two different ways below.

Indeed, a simple measurement in the computational basis will only provide us with estimations of the probabilities that are the squares of the quantum amplitudes (see Section A.1). In the case of neural networks, it is important to obtain the sign of the layer’s components in order to apply certain type of non linearities. For instance, the ReLu activation function is often used to set all negative components to 0.

In Fig.9, we propose a specific enhancement to our circuit to obtain the signs of the vector’s components at low cost. The sign retrieval procedure consists of three parts.

1. a)

   The circuit is first applied as described above, allowing to retrieve each squared amplitude $y_{j}^{2}$ with precision $\delta>0$ using the $\ell_{\infty}$ tomography. The probability of measuring the unary state $\mathinner{|{e_{1}}\rangle}$ (i.e. $\mathinner{|{100...}\rangle}$), is $p(e_{1})=y_{1}^{2}$.

2. b)

   We apply the same steps a second time on a modified circuit. It has additional RBS gates with angle $\pi/4$ at the end, which will mix the amplitudes pair by pair. The probabilities to measure $\mathinner{|{e_{1}}\rangle}$ and $\mathinner{|{e_{2}}\rangle}$ are now given by $p(e_{1})=(y_{1}+y_{2})^{2}$ and $p(e_{2})=(y_{1}-y_{2})^{2}$. Therefore if $p(e_{1})>p(e_{2})$, we have $sign(y_{1})\neq sign(y_{2})$, and if $p(e_{1})<p(e_{2})$, we have $sign(y_{1})=sign(y_{2})$. The same holds for the pairs ($y_{3}$, $y_{4}$), and so on.

3. c)

   We finally perform the same where the RBS are shifted by one position below. Then we compare the signs of the pairs ($y_{2}$, $y_{3}$), ($y_{4}$, $y_{5}$) and so on.

At the end, we are able to recover each value $y_{j}$ with its sign, assuming that $y_{1}>0$ for instance. This procedure has the benefit of not adding depth to the original circuit, but requires 3 times more runs. The overall cost of the tomography procedure with sign retrieval is given by $\widetilde{O}(n/\delta^{2})$.

In Fig.10 we propose another method to obtain the values of the amplitudes and their signs, which is in fact what we used for the hardware demonstrations. Compared to the above procedure, it relies on one circuit only, but requires an extra qubit and a depth of $3n+O(1)$ instead of $2n+O(1)$.

This circuit performs a Hadamard and CNOT gate in order to initialize the qubits in the state $\frac{1}{\sqrt{2}}\mathinner{|{0}\rangle}\mathinner{|{0}\rangle}+\frac{1}{\sqrt{2}}\mathinner{|{1}\rangle}\mathinner{|{e_{1}}\rangle}$, where the second register corresponds to the $n$ qubits that will be processed by the pyramidal circuit and the loaders.

Next, applying the data loader for the normalized input vector $x$ (see Section 2.3) and the pyramidal circuit will, according to Eq.(4), map the state to

In other words, we performed the pyramid circuits controlled on the first qubit being in state $\mathinner{|{1}\rangle}$. Then, we flip the fisrt qubit with an X gate and perform a controlled loading of the uniform norm-1 vector $(\frac{1}{\sqrt{n}},\cdots,\frac{1}{\sqrt{n}})$. For this, we add the adjoint data loader for the state, a CNOT gate and the data loader a second time. Recall that if a circuit $U$ is followed by $U^{\dagger}$, it is equivalent to the identity. Therefore, this will load the uniform state only when the first qubit is in state $\mathinner{|{1}\rangle}$:

Finally, a Hadamard gate will mix both parts of the amplitudes on the extra qubit to give us the desired state:

On this final state, we can see that the difference in the probabilities of measuring the extra qubit in state $0$ or $1$ and rest in the unary state $e_{j}$ is given by $\Pr[0,e_{j}]-\Pr[1,e_{j}]=\frac{1}{4}\left(W_{j}x+\frac{1}{\sqrt{n}}\right)^{2}-\frac{1}{4}\left(W_{j}x-\frac{1}{\sqrt{n}}\right)^{2}=W_{j}x/\sqrt{n}$. Therefore, for each $j$, we can deduce the sign of $W_{j}x$ by looking at the most frequent output of the measurement of the first qubit. To deduce as well the value of $W_{j}x$, we simply use $\Pr[0,e_{j}]$ or $\Pr[1,e_{j}]$ depending on the sign found before. For instance, if $W_{j}x>0$ we have $W_{j}x=2\sqrt{\Pr[0,e_{j}]}-\frac{1}{\sqrt{n}}$.

Combining with the $\ell_{\infty}$ tomography and the non linearity, the overall cost of this tomography is given by $\widetilde{O}(n/\delta^{2})$ as well.

In the previous sections, we have seen how to implement a quantum circuit to perform the evolution of one orthogonal layer. In classical deep learning, such layers are stacked to gain in expressivity and accuracy. Between each layer, a non-linear function is applied to the resulting vector. The presence of these non-linearities is key in the ability of the neural network to learn any function [15].

The benefit of using our quantum pyramidal circuit is the ability to simply concatenate them to mimic a multi layer neural network. After each layer, a tomography of the output state $\mathinner{|{z}\rangle}$ is performed to retrieve each component, corresponding to its quantum amplitudes (see Section 3.2). A non linear function $\sigma$ is then applied classically to obtain $a=\sigma(z)$. The next layer starts with a new unary data loader (See Section 2.3). This hybrid scheme allows as well to keep the depth of the quantum circuits reasonable for NISQ devices, by applying the neural network layer by layer.

Note that the quantum neural networks we propose here are close to the behaviour of classical neural networks and thus we can control and understand the quantum mapping and implement each layer and its non linearities in a modular way. They are also different regarding the training strategies which are close to the classical ones but utilise a different optimization landscape that can provide diffrent models (see Section 4.2 for details). It will be interesting to compare our pyramidal circuit to a quantum variational circuit with $n$ qubits and $n(n-1)/2$ gates of any type, as we usually see in the literature. Using such circuits we would explore among all possible $2^{n}\times 2^{n}$ matrices instead of $n\times n$ classical orthogonal matrices, but so far there’s no theoretical ground to explain why this should provide an advantage.

As an open outlook, one could imagine incorporating additional entangling gates after each pyramid layer (composed, for instance, of $CNOT$ or $CZ$). This would mark a step out of the unary basis and could effectively allow exploring more interactions in the Hilbert Space.

The backpropagation in a fully connected neural network is a well known and efficient procedure to update the weight matrix at each layer [16, 17]. At layer $\ell$, we note its weight matrices $W^{\ell}$ and biases $b^{\ell}$. Each layer is followed by a non linear function $\sigma$, and can therefore be written as

After the last layer, one can define a cost function $\mathcal{C}$ that compares the output to the ground truth. The goal is to calculate the gradient of $\mathcal{C}$ with respect to each weight and bias, namely $\frac{\partial\mathcal{C}}{\partial W^{\ell}}$ and $\frac{\partial\mathcal{C}}{\partial b^{\ell}}$. In the backpropagation, we start by calculating these gradients for the last layer, then propagate back to the first layer.

We will require to obtain the error vector at layer $\ell$ defined by $\Delta^{\ell}=\frac{\partial\mathcal{C}}{\partial z^{\ell}}$. One can show the backward recursive relation $\Delta^{\ell}=(W^{\ell+1})^{T}\cdot\Delta^{\ell+1}\odot\sigma^{\prime}(z^{\ell})$, where $\odot$ symbolizes the Hadamard product, or entry-wise multiplication. Note that the previous computation requires simply to apply the layer (ie apply matrix multiplication) in reverse. We can then show that each element of the weight gradient matrix at layer $\ell$ is given by $\frac{\partial\mathcal{C}}{\partial W^{\ell}_{jk}}=\Delta^{\ell}_{j}\cdot a^{\ell-1}_{1}$. Similarly, the gradient with respect to the biases is easily defined as $\frac{\partial\mathcal{C}}{\partial b^{\ell}_{j}}=\Delta^{\ell}_{j}$.

Once these gradients are computed, we update the parameters using the gradient descent rule, with learning rate $\lambda$:

Looking through the prism of our pyramidal quantum circuit, the parameters to update are no longer the individual elements of the weight matrices directly, but the angles of the RBS gates that give rise to these matrices. Thus, we need to adapt the backpropagation method to our setting based on the angles. We will start by introducing some notation for a single layer $\ell$, which will not be explicit in the notation for simplicity. We assume we have as many output bits as input bits, but this can easily be extended to the rectangular case.

We first introduce the notion of timesteps inside each layer, which correspond to the computational steps in the pyramidal structure of the circuit (see Fig.13). It is easy to show that for $n$ inputs, there will be $2n-3$ such timesteps, each one indexed by an integer $\lambda\in[0,\cdots,\lambda_{max}]$. Applying a timestep consists in applying the matrix $w^{\lambda}$, made of all the RBS gates aligned vertically at this timestep ($w^{\lambda}$ is the unitary in the unary basis, see Section 3 for details). Each time a timestep is applied, the resulting state is a vector in the unary basis named inner layer and noted $\zeta^{\lambda}$. This evolution can be written as $\zeta^{\lambda+1}=w^{\lambda}\cdot\zeta^{\lambda}$. We use this notation similar to the real layer $\ell$, with the weight matrix $W^{\ell}$ and the resulting vector $z^{\ell}$ (see Section 4.1).

In fact we have the correspondences $\zeta^{0}=a^{\ell-1}$ for the first inner layer, which is the input of the actual layer, and $z^{\ell}=w^{\lambda_{max}}\cdot\zeta^{\lambda_{max}}$ for the last one. We also have $W^{\ell}=w^{\lambda_{max}}\cdots w^{1}w^{0}$.

We use the same kind of notation for the backpropagation errors. At each timestep $\lambda$ we define an inner error $\delta^{\lambda}=\frac{\partial\mathcal{C}}{\partial\zeta^{\lambda}}$. This definition is similar to the layer error $\Delta^{\ell}=\frac{\partial\mathcal{C}}{\partial z^{\ell}}$. In fact we will use the same backpropagation formulas, without non linearities, to retrieve each inner error vector $\delta^{\lambda}=(w^{\lambda})^{T}\cdot\delta^{\lambda+1}$. In particular, for the last timestep, the first to be calculated, we have $\delta^{\lambda_{max}}=(w^{\lambda_{max}})^{T}\cdot\Delta^{\ell}$. Finally, we can retrieve the error at the previous layer $\ell-1$ using the correspondence $\Delta^{\ell-1}=\delta^{0}\odot\sigma^{\prime}(z^{\ell})$.

The reason for this breakdown into timesteps is the ability to efficiently obtain the gradient with respect to each angle. Let’s consider the timestep $\lambda$ and one of its gate with angle noted $\theta_{i}$ acting on qubits $i$ and $i+1$ (note that the numbering is different from Fig.13). We will decompose the gradient $\frac{\partial\mathcal{C}}{\partial\theta_{i}}$ using each component, indexed by the integer $k$, of the inner layer and inner error vectors $\frac{\partial\mathcal{C}}{\partial\theta_{i}}=\sum_{k}\frac{\partial\mathcal{C}}{\partial\zeta^{\lambda+1}_{k}}\frac{\partial\zeta^{\lambda+1}_{k}}{\partial\theta_{i}}=\sum_{k}\delta^{\lambda+1}_{k}\frac{\partial(w^{\lambda}_{k}\cdot\zeta^{\lambda})}{\partial\theta_{i}}$, where $w^{\lambda}_{k}$ is the k^th row of matrix $w^{\lambda}$.

Since timestep $\lambda$ is only composed of separated RBS gates, the matrix $w^{\lambda}$ consists of diagonally arranged $2\times 2$ block submatrices given in Eq.(1). Only one of these submatrices depends on the angle $\theta_{i}$ considered here, at the position $i$ and $i+1$ in the matrix. We can thus rewrite the above gradient as $\frac{\partial\mathcal{C}}{\partial\theta_{i}}=\delta^{\lambda+1}_{i}\frac{\partial}{\partial\theta_{i}}\left(w^{\lambda}_{i}\cdot\zeta^{\lambda}\right)+\delta^{\lambda+1}_{i+1}\frac{\partial}{\partial\theta_{i}}\left(w^{\lambda}_{i+1}\cdot\zeta^{\lambda}\right)$, or:

Therefore we have shown a way to compute each angle gradient: During the feedforward pass, one must apply sequentially each of the $2n-3=O(n)$ timesteps, and store the resulting vectors, the inner layers $\zeta^{\lambda}$. During the backpropagation, one obtains the inner errors $\delta^{\lambda}$ by applying the timesteps in reverse. One can finally use a gradient descent on each angle $\theta_{i}$, while preserving the orthogonality of the overall equivalent weight matrix $\theta_{i}^{\ell}\leftarrow\theta_{i}^{\ell}-\lambda\frac{\partial\mathcal{C}}{\partial\theta_{i}^{\ell}}$. Since the optimization is performed in the angle landscape, and not on the equivalent weight landscape, it can potentially be different and produce different models. We leave open the study of the properties of both landscapes.

As one can see from the above description, this is in fact a classical algorithm to obtain the angle’s gradients, which allows us to train our OrthoNN efficiently classically while preserving the strict orthogonality. To obtain the angle’s gradient, one needs to store the $2n-3$ inner layers $\zeta^{\lambda}$ during the feedforward pass. Next, given the error at the following layer, we perform a backward loop on each timestep (see Fig.12). At each timestep, we obtain the gradient for each angle parameter, by simply applying Eq.(12). This requires $O(1)$ operations for each angle. Since there are at most $n/2$ angles per timesteps, estimating gradients has a complexity of $O(n^{2})$. After each timestep, the next inner error $\delta^{\lambda-1}$ is computed as well, using at most $4n/2$ operations.

In the end, our classical algorithm allows us to compute the gradients of the $n(n-1)/2$ angles in time $O(n^{2})$, thus performing a gradient descent respecting the strict orthogonality of the weight matrix in the same time. This is considerably faster than previous methods based on Singular Value Decomposition methods and provides a training method that is asymptotically as fast as for normal neural networks, while providing the extra property of orthogonality.