We present, , a novel layer-wise stochastic optimiser tailored towards training neural networks with binary weights, known as binary neural networks (BNNs), on quantum hardware. BNNs reduce the computational requirements and energy consumption of deep learning models
with minimal loss in accuracy. However, training them in practice remains to be an open challenge. Most known BNN-optimisers either rely on projected updates or binarise weights post-training.
Instead,  approximately maps the gradient onto binary variables, by solving a quadratic constrained binary optimisation. Under practically reasonable assumptions, we show that this update rule converges with a rate of $\mathcal{O}(1 / \sqrt{T})$. Moreover, we show how the $\mathcal{NP}$-hard projection can be effectively executed on an adiabatic quantum annealer, harnessing recent advancements in quantum computation. We also introduce a projected version of this update rule and prove that if a fixed point exists in the binary variable space, the modified updates will converge to it. Last but not least, our algorithm is implemented layer-wise, making it suitable to train larger networks on resource-limited quantum hardware. Through extensive evaluations, we show that  outperforms or is on par with competitive and well-established baselines such as BinaryConnect, signSGD and ProxQuant when optimising the Rosenbrock function, training BNNs as well as binary graph neural networks.

§ INTRODUCTION

Our contemporary times witness the emergence of two exciting and prominent avenues of scientific advancement: quantum computing (QC) [26, 51] in physics and
machine learning [40, 80] in computer science. QC seeks to develop novel, efficient computing systems that can solve problems beyond the capabilities of classical computers. On the other hand, machine learning aims to create algorithms that can analyze and learn from data without the need for explicit programming. With these advancements in mind, it is natural to pose the question of whether quantum computers can be utilised to train our learning machines or in other words, to optimise their parameters.

This seemingly simple question calls for a
reconsideration of our training algorithms, due to an important distinction between classical optimisation and the state-of-the-art quantum annealers (QA), such as the D-Wave Advantage [22, 12][Quantum annealers based on the adiabatic quantum computation have demonstrated remarkable progress in the scale of quantum computing, reaching beyond 5000 qubits, compared to the universal gate machines such as Google Sycamore (53 qubits) and IBM Eagle (127 qubits).]. While classical optimisation algorithms commonly use relax-and-round schemes to tackle discrete problems, QAs excel in solving (potentially non-convex) combinatorial, binary optimisation problems that can be expressed in the form of Ising models [45]. Therefore, any training algorithm looking to leverage QA should be modified to involve (sub-)problems in the form of an Ising model. For standard training algorithms, this is not the case.

In this work, we fill this gap for the case of binary (graph) neural networks (BGNNs) [54, 68].
Due to their binary nature, training BGNNs effectively involves solving a combinatorial optimisation problem, which is supported by real quantum hardware such as quantum annealers (QAs). Once trained, BGNNs can be deployed on edge devices for efficient inference.

To train BGNNs on quantum hardware, we introduce , a hybrid classical-quantum, stochastic projected-gradient descent optimiser where the gradients are approximated by their projections onto the set of binary variables by solving a Quadratic Unconstrained Binary Optimisation (QUBO) problem [45] on a real quantum annealer. In contrast to projection-based optimisers [54, 6],  respects the binary nature of the parameter space and potentially delivers higher quality weight updates. Our generalised algorithm is developed in two stages. We develop (i) , a continuous optimiser where the gradients are replaced by their binary projections, and (ii) , a projected, binary version where the parameters are updated taking into account their binary nature at all times.
Drawing connections to the state-of-the-art binary gradient optimisers such as signSGD [58], we ensure that the convergence rate of (i) is $\mathcal{O}(1 / \sqrt{T})$, and for (ii), we show that if a fixed point exists in the space of binary variables, our algorithm will converge to it. We empirically validate our assumptions required to prove these results, which we demonstrate are milder compared to the assumptions made in existing algorithms such as signSGD [58]. This remarkably justifies why rethinking our optimisers in the context of quantum hardware is not a mere implementation concern but a way to make progress in the accuracy of training algorithms.

The most computationally expensive step in our algorithm is the projection onto the non-convex set of binary variables, which adheres to a QUBO form and hence can be computed on a quantum annealer. As quantum hardware resources are limited, closed-form, one-shot approaches such as that of Sasdelli  [59] fail to train even modest-sized neural networks. As a remedy, we show how to deploy  scalably, by computing weight updates layer-wise and incrementally over data batches, inspired by the recent layer-wise optimisers [74]. This is illustrated in <ref>. Finally, we extend our formulation to binary graph-convolutional neural networks [39, 68] and conduct extensive evaluations on both image and graph representation learning.
In summary, our contributions are:

* We propose , a novel, stochastic optimiser tailored for training binary (graph) neural networks utilizing real quantum hardware.

* We prove, under mild assumptions, that the continuous version of our algorithm, , converges with a rate of $\mathcal{O}(1 / \sqrt{T})$ using binary-projected gradients.

* We prove that additionally projecting the weights onto binary variables still achieves convergence as long as a fixed point exists in the set of binary variables.

* We show an equivalence of our binary projection to a specific QUBO problem, allowing us to implement our algorithm on quantum hardware.

* We conduct thorough evaluations for all variants of our algorithm on the actual quantum annealer of D-Wave as well as on classical simulated annealers. These experiments reveal that our method (i) can outperform classical benchmarks in both settings and (ii) can train binary neural networks in a scalable manner on quantum annealers.

To the best of our knowledge, this is the first time a general, practical binary (graph) neural network is trained on an actual adiabatic quantum computer. Our algorithm is intended to harness the current and upcoming advancements in QC technology, effectively circumventing practical hardware limitations such as those imposed by memory constraints. We will make our implementation available upon publication.

§ RELATED WORK

Training binary neural networks
Compared to full-precision neural networks, BNNs [33] consume less memory and provide faster inference at the cost of certain information loss. Much effort has been devoted to mitigating performance degradation due to binarisation, as detailed in recent surveys [55, 75].
Early, naive methods such as BinaryConnect (BC) [21], BinaryNet [33], BinaryNetg [33] or Dorefa-Net [81] aim to binarise weights, activations and during inference and training. Optimisation-based BNNs such as XNOR-Net [56], XNOR-Net++ [15], Parameterised Clipping Activation (PACT) [18], LAB [31], HWGQ-Net [16] or ReActNet [43] attempt to more directly address the accuracy drop resulting from weight and activation binarisation in the naive models, while preserving the compact nature [55].
Later, ProxQuant [6] approached binarisation by using a gradual regulariser, which starts with continuous weights and during training gradually regularises those to binary ones.
Standard techniques are used in optimisation, involving Adam [38], AdaMax [38], SGD and RMSprop [24]. We refer the reader to Qin et al. [55] for a survey.
The approach of [5] provides a binarisation scheme for GNNs relying on XNOR-Net++ and knowledge transfer with the student-teacher model.
Because quantum annealers operate on binary variables, instead of using an approximation/binarisation scheme, we directly calculate the binary weight updates for the networks.

Quantum neural networks
As an emerging field, quantum machine learning [10, 61, 14, 71] has shown that the training of linear regression, support vector machines and k-means clustering admit QUBO-like formulations [3]. The first neural network variants trained on quantum hardware were Boltzmann Machines (BMs) as they naturally lent themselves to quantum annealing [27, 10, 70, 1]. Quantum deep learning [70, 29, 36], i.e., creating quantum circuits that enhance the operations of neural networks by physically realizing them, has emerged to alleviate some of the computational limitations of classical deep learning, thanks to the efficient training algorithms [37, 7]. However, a severe challenge remains to be implementing non-linearities and other non-unitary operations with quantum unitaries [60, 36, 57]. Similarly, quantum convolutional neural networks (Q-CNNs) [20] and their graph counterparts (QGCNs) [66] have been developed with the mindset of implementing the analogue of classical CNN/GNN operations via quantum gates. This is different from our approach where we leverage QC as a computational tool to train classical architectures.
Despite the widespread use of quantum annealers in various related domains such as quantum computer vision [62, 11, 76], training neural networks via QAs is less attended. Only recently, Sasdelli  [59] proposed a one-shot, closed-form way to train binary neural networks but this approach is not practically applicable to even small-sized neural networks. We circumvent this problem by adopting a layerwise approach, updating each layer individually, leading to a more scalable training algorithm.

§ METHOD

We consider a differentiable, Lipschitz continuous function $E_f(\x) := E(f(\x)) : \R^n \rightarrow \R$ defined on $f: \R^n \rightarrow \R^m$ and the loss function $E:\R^m \rightarrow \R$, where $\y = f(\x)$ denotes the evaluation of $f$ at $\x$ and $n$ is the dimension of the input of $f$. Let $K_i$ denote the Lipschitz constant for the $i^{th}$ dimension of $E_f$ and $\bar{K} = \frac{1}{n} \sum_{i = 1}^nK_i$. While $t$ denotes the current iteration, $T$ is the number of total iterations. We define $\{\alpha_t\}_t$ as a series of learning rates, $\{ \x^t\}_t$ as the series of iterates and $\tilde{\nabla}_{\x}E_f(\x)$ refers to the stochastic gradient w.r.t. $\x$.
Further let $\bZ^t\in\R^{n\times m}$ whose $i^\mathrm{th}$ column is $\R^n\ni\bz_i^t = \bj^t_i / \| \bj_i^t\|^2 $ where $\bj_i^t = \nabla_{\x} y_i|_{\x^t}$. Note that, $\bj_i^t$ denotes the $i^\mathrm{th}$ column of the Jacobian map $\R^{n\times m}\ni\J^t:= \nabla_{\x} \y$ evaluated at $\x^t$.
We further define $\bU=[\bu_i]_i$ and let $\bv \in \R^m$.
Our goal is to design an algorithm for training binary (graph) neural networks which is:

* QA-friendly: utilizing quantum hardware to update the weights while ensuring that the results obtained are accurate and reliable, through a hybrid feedback loop;

* Provably convergent: training should result in the network reaching a stable and optimal state;

* Layer-wise compliant: being able to train individual layers of the network independently, thereby enabling scalable training of large and complex models.

To meet these criteria, we introduce a new optimiser,  (), that maintains the binary nature of weights at all times by projecting the gradients (weight updates) onto the set of binary variables.
Our approach begins with  (), a continuous optimiser that approximates real-valued gradients with binary ones. Under mild assumptions, which we verify in our supplementary material, this optimiser can converge to a local, fixed point[In our context, a fixed point differs from a local minimum. A given binary gradient (scaled with the learning rate) may not be sufficient for jumping out of a fixed point regardless of the Hessian.]. We then introduce , a projected version of  that can converge to a fixed point in the binary variable space if such a fixed point exists. Finally, we will show how to express our projection as a Quadratic Unconstrained Binary Optimisation (QUBO), solvable on an adiabatic quantum computer (such as D-Wave) resulting in Quantum   (Q-).
In what follows, we will provide proof sketches and leave the full proofs to our supplementary material.

§.§  ()

We start by presenting our new binary mapping operator before defining .
[Binary map (BM)] We define
$\bPi_{\bU}:\R^m \rightarrow \{ \pm 1\}^{n}$ to be the map

\begin{align}\label{eq:BP_paper}
\bPi_{\bU} (\bv) := \argmin_{\bg \in \{-1,1 \}^{n}} \sum_{i=1}^m||v_i - \bg^\top \bu_i||_2^2.
\end{align}

We also let $\hat{\bPi}_{\bU} :\R^m\to\R^n$ denote the relaxed or continuous version of our projection map. A direct application of the chain rule reveals that such a map satisfies:

\begin{align}
\label{eq:binary_relax}
\hat{\bPi}_{\Z^t}(\tilde{\nabla}_{\y}E_f(\x)) = \argmin_{\bb\in\R^{n}} \|\tilde{\nabla}_{\x}E_f(\x) |_{\x^t} -\bb\|_2^2.
\end{align}

Intuitively, $\hat{\bPi}$ acts as a Jacobian, transforming the tangent plane of the gradient with respect to $\y$ onto the tangent plane of the gradient with respect to $\x$.
However, our original operator $\bPi$ projects onto the binary numbers and not the reals. This non-convex map only approximates $\hat{\bPi}$. Hence, we write:

\begin{align}
\bPi_{\bZ^t}(\tilde{\nabla}_{\y}E_f |_{\x^t}) \approx \argmin_{\bb\in\{-1,1\}^{n}} \|\tilde{\nabla}_{\x}E_f |_{\x^t} -\bb\|_2^2.
\end{align}

Our supplementary material shows that our convergence depends not on the quality of this approximation but rather on the direction in which the binary projection and the real gradient are pointing.

Note that, usually $m\leq n$ rendering <ref> under-determined. Yet, in practice, as we operate on the confined space of binary variables we get valid solutions, one of which can approximate the gradient w.r.t. the input.

With this motivation, we use $\bPi_{\Z^t}(\tilde{\nabla}_{\y}E_f(\x))$ as an approximation for the true stochastic gradient with respect to $\x$ and introduce our new optimiser,  ():

The  algorithm admits the following update rule:

\begin{align}
\x^{t + 1} = \x^t - \alpha_t\bPi_{\Z^t}(\tilde{\nabla}_\by E_f(\bx^t)).
\end{align}

This update rule resembles signSGD [58] up to a difference in projection.
We are able to show that  will converge to a stationary optimum under the following conditions:
[Lower bound]
There exits a value $E_f^* \in \R$ with $E_f^* \leq E_f(\x), \quad \forall \x \in \R^n$.

[Co-Directionality Probability (CDP)]
The probability that the binary gradient approximation of BM points in the same direction as the true gradient at the point $\x^t$ is larger than $1/2$:

\begin{align}
\label{eq:defSPB}\tag{CDP}
&= \operatorname{Prob}\left( \bPi_{\Z^t}(\tilde{\nabla}_\y E_f(\x^t))_i=\operatorname{sign} (\nabla_\x E_f(\x^t))_i\right)>\frac{1}{2}. \nonumber
\end{align}

CDP is similar to the success probability (SP) assumption in signSGD [58], yet the empirical evidence we provide in our supplementary material demonstrates that this condition is most likely satisfied and that the SP of our binarisation is higher than the one of signSGD.

To measure our gradient in relation to the CDP, we use $\rho$-norm from Safaryan and Richtárik [58].
[Gradient $\rho$-norm [58]]

Let $\rho_i$ be the probability function from the CDP assumption. We define the $\rho$-norm of gradient $\bg$ as follows:

\begin{align}
\|\bg(\x)\|_\rho:=\sum_{i=1}^n\left(2 \rho_i(\x)-1\right)\left|g_i(\x)\right|.
\end{align}

If CDP holds, then <ref> is a norm in the mathematical sense,  fulfills subadditivity, homogeneity and positive definiteness.

Our proof relies on the proof of [58], where we replace SP by CDP.

We are now ready to state the convergence of .

Assuming <ref>, <ref>, and using a step size of $\alpha_t=\alpha_0 / \sqrt{t+1}$,  guarantees a convergence rate of $\mathcal{O}({1}/{\sqrt{T}})$:

\begin{align}
\begin{split}
&\min _{0 \leq t<T} \mathbb{E}\left\|\nabla_{\x} E_f\left(\x^t\right)\right\|_\rho \\
&\leq \frac{E_f\left(\x^0\right)-E_f^*}{\alpha_0 \sqrt{T}}+\frac{3 \alpha_0 n \bar{K}}{2} \frac{\log T}{\sqrt{T}}.
\end{split}
\end{align}

If $\alpha_t \equiv \alpha>0$, we get a convergence of $\mathcal{O} ({1}/{T})$ convergence to a neighbourhood:

\begin{align}
\begin{split}
\frac{1}{T} \sum_{t=0}^{T-1} \mathbb{E}\left\|\nabla_{\x} E_f\left(\x^t\right)\right\|_\rho \leq \frac{E_f\left(\x^0\right)-E_f^*}{\alpha T}+\frac{\alpha n \bar{K}}{2}.
\end{split}
\end{align}

The proof follows very closely the proof of the convergence of signSGD  [58] as well as the standard SGD proof. The only difference lies in the expected value of the interaction of our projected stochastic gradient and the real gradient. We, therefore, show that the following holds:

\begin{align*}
\mathbb{E}\left[\left\langle \nabla_{\bx}E_f(\x^t), \bPi_{\Z^t}( \widetilde{\nabla}_{\y}E_f(\x^t))\right\rangle\right]=\left\|\nabla_{\bx}E_f(\x^t)\right\|_\rho. \nonumber\vspace{3cm}\tag*{\qedhere}
\end{align*}

§.§  ()
The scheme we introduced in <ref> updates the parameters of a function in the real values.
We, however, ultimately introduce an optimiser for training binary (graph) neural networks.
Hence, we post-project the updates and introduce :

We now devise a projected variant of SBGD with the distinction that we evaluate the gradients on the variables restricted to $\{ \pm 1\}^n$:

\begin{align}\label{eq:BAPGDupdate}
\begin{split}
\hat{\x}^t &=\sign(\x^t) \\
\x^{t+1}&=\x^t-\alpha_t \bPi_{\Z^t}(\widetilde{\nabla}_{\y} E_f\left(\hat{\x}^t\right)).
\end{split}
\end{align}

This approach bears similarities to a variety of algorithms in the literature. First, <ref> and resembles ProxQuant [6] with the critical distinction that while ProxQuant employs the continuous gradient, our approach utilises the binarised gradient directly. Second, by drawing connections to manifold optimisation [13], our projection parallels a Euclidean-to-Riemannian gradient operator whereas the projected update mirrors an approximate retraction. Yet, this is merely an analogy as we operate in a discrete space.

We now show that if a fixed point exists in the binary set of variables that minimises our objective,  algorithm can converge to such a point.
To this end, we will further assume that the sum of our learning rates will be infinite:
[Divergent learning rates]
The learning rates $\alpha_t$ do not converge, hence $\sum_t \alpha_t = \infty$. Otherwise, we directly could assume the convergence of the algorithm.

Under <ref>, $\mathbf{s} \in \{\pm 1\}^n$ is a fixed point for <ref> if and only if $\sign (\bPi_{\mathbf{s}} (\nabla_{\mathbf{y}}E_f(\mathbf{s}))_i) = - \mathbf{s}_i$. This point might not exist, in which case  does not converge.

The proof follows analogously to the proof of Prop. 5.3 in Bai et al. [6].

§.§  ()

Finding an ideal binary solution for the optimisation problem stated in <ref> is an $\mathcal{NP}$-hard problem.
Therefore, we opt for a Quantum Annealer to effectively solve this problem.
To deploy the binary-projection problem onto a QA, we construct the corresponding Quadratic Unconstrained Binary Optimisation (QUBO) problem:

The binary projection $\bPi_{\Z^t}(\bv)$ in <ref> admits the following Ising Model or quadratic unconstrained binary optimisation (QUBO) form:

\begin{align}\label{eq:BP_as_qubo_paper}
\bPi_{\bU} (\bv) &= \argmin_{\bg \in \{-1,1 \}^{n}} \bg^\top\sum_{i=1}^m\bQ_i\bg + \bs^\top\bg
\end{align}

\begin{align}
\bs=-2\sum\limits_{i=1}^m v_i\bu_i^t,\quad \bQ_i={\bu_i^t}{\bu_i^t}^\top \quad\mathrm{and}\quad \bQ=\sum\limits_{i=1}^m\bQ_i.\nonumber
\end{align}

The derivation follows from the expansion of the Euclidean norm.

<ref> allows us to compute the costly projection step on a quantum annealer such as D-Wave [23, 22].
To compute the projection for a batch of data, we are using the multi-objective optimisation of [4]. With that, we define the quantum versions of our optimisers as follows:
[, ]
Whenever we compute the projection by solving the QUBO given in <ref>, we will use the prefix Q-, to denote the quantum-implementable variants of  and .

§.§ Layer-wise Training

Training data $\cD = \{ (\bx_i, \hat{y}_i) \}^{D}_{i=1}$, batch size $B$, learning rate $\alpha$, real initial weights $\{ {\mathbf{\Omega}}^{\ell} \}_{\ell = 0}^{L-1} $
$t \in [1, \dots, T]$
$\{ \bW^{\ell} \}_{\ell = 1}^{L} \leftarrow \sign(\bOmega^\ell)$
Sample a batch index set $\cB \subset \{1, \dots, D\}$.
$ \by_{\cB} \leftarrow$ Feedforward pass of $\x_\cB$ with <ref>.
$\{ \dot{\br}^{\ell}_{i,\cB} \}_{\ell = 1}^L \leftarrow$ Compute
intermediate gradients
2.5em by <ref> for training data
$\ell = 1,\dots,L$
${\dot{\bW}}^\ell \leftarrow [\bPi_{\Z^{t,\ell}_{\mathcal{B}, i}}(\dot{\br}^\ell_{i,\mathcal{B}})]_{i = 1}^m$ By solving the QUBO  4.5em defined  in <ref>
${\mathbf{\Omega}}^{\ell} \leftarrow {\mathbf{\Omega}}^{\ell} - \alpha {\dot{\bW}}^\ell$

We now apply  layer-wise and to the training of binary graph-convolutional neural networks (GCNs) [39]. As MLPs are a special case of binary GCNs this method, summarised in <ref>, can also be used to train binary MLPs.
We start by considering a general binary neural network where $\R^{m\times n}\ni\X:=\{\x_i\}_{i\leq n}$ denotes the input features and $\X^\ell$ to be the input of the $\ell^{\text{th}}$ layer.
We denote the real-valued weight matrices as
$\{ \mathbf{\Omega}^\ell\}_{\ell=1}^L$ where $\mathbf{\Omega}^\ell \in \R^{n_{\ell}\times m_{\ell}}$. Similarly, $\W^{\ell} := \sign \mathbf{\Omega}^{\ell}$ denote layer-wise binary weights.
We first show how to update the weights of each layer, individually.
We further update each column $j$ of the weight matrix individually.
Let therefore be $\bZ^{t,\ell}_{\mathcal{B},j}\in\R^{n\times m}$ whose $i^\mathrm{th}$ column is $\R^n\ni\bz_i^t = \bj^t_i / \| \bj_i^t\|^2 $ where $\bj_i^t = \frac{\partial \br_i}{\partial \bw_j}|_{\X_{\mathcal{B}}}$. Note that, $\bj_i^t$ denotes the $i^\mathrm{th}$ column of the Jacobian map $\R^{n\times m}\ni\J^t:= \frac{\partial \y}{\partial \x}$ evaluated at $\X_{\mathcal{B}}$ with the weights $\W^\ell$.

[Intermediate gradient]

Let $\cB \subset [1\dots D]$ denote the index for a batch of input vectors. We denote the output vectors of the $\ell^{\mathrm{th}}$ layer corresponding to the batch $\cB$ as $\br_\cB^{\ell}$ and define:

\begin{align}
\bR^{\ell}_\cB = ((\mathsf{h} \cdots \; \mathsf{h}({\X}\W^{1}) \cdots {\W}^{\ell-1} ){\mathbf{\Omega}}^{\ell})_\cB.
\label{eq:pseudo_label}
\end{align}

where the binary weights $\W^{\ell} := \sign \mathbf{\Omega}^{\ell}$ pass the information through the non-linear activation function $\mathsf{h}(\X) =\mathsf{hardTanh}(\X)$ [19].
This is the straight-through estimator (STE), also used in [64, 8].
The intermediate gradient $\dot{\bR}$ is then defined as:
$\dot{\bR}^{\ell}_\cB = \frac{\partial E}{\partial \bR^{\ell}_\cB}$, whose $i^{th}$ row is $\br^\ell_{i,\cB}$.

We now extend this definition for a general binary MLP to the case of graph-convolutional networks.
[Binary GCN (b-GCN)]
We consider an input graph $\Graph=(\Vertex=\{1,2,\dots,n\}, \Edge\subset \Vertex\times\Vertex)$ of $n$ vertices. We override $\R^{m\times n}\ni\X:=\{\x_i\}_{i\leq n}$ to denote node features associated to each and every vertex. We encode the edge structure of the graph in a adjacency matrix $\A\in\R^{n\times n}$ and let $\hat{\A}=\A+\Id$.
The forward pass is defined as:

\begin{align}
\small
\y_j \label{eq:forward_pass} = \mathrm{softmax}\left(\hat{\A}\mathsf{h}\left( \cdots \;
\mathsf{h}\left({\hat{\A}\X}\W^{1}\right) \cdots \right)
{\mathbf{W}}^{L}\right)_j. %
\end{align}

It is common in GCNs to feed the connected components as a whole through the network to ensure the correctness of convolutions.

§ EXPERIMENTS

In this section, we investigate how well  performs on various problem types.
We start by showing the minimisation of the Rosenbrock curve.
Next, we assess the convergence in the binary logistic regression problem setup and evaluate the accuracy as well as convergence of  on binary MLPs.
Subsequently, we train and test binary graph-convolutional networks (GCNs) on two different datasets.
Finally, we study the Hamiltonian gap of our QUBO problems to observe their difficulty.
For our method we are either using the real quantum annealer developed by D-Wave or the well-engineered simulated annealing software, Gurobi. Using D-Wave's QA is marked with (D) and Gurobi's SA with (G).
In the binary setups (all except the Rosenbrock function in <ref>), we compare our method against the strong baselines of BinaryConnect using either Stochastic Gradient Descent (BC (SGD))  [21] or signSGD (BC (signSGD)) [21, 58] as well as against ProxQuant [6] using the SGD optimiser without momentum.
In all figures, the shading captures the variance over five runs and the learning rates have been tuned for each optimiser.
We provide additional details, ablation studies and training of a continuous neural network in our supplementary material.

§.§ Minimisation of the Rosenbrock function

Minimisation of the Rosenbrock function; $y$-axis is in $\log$-scale.

We start by assessing  in the classical setting of finding the minimum of the Rosenbrock function, $f(\x)=100\left(x_{2}-x_1^2\right)^2+\left(1-x_1\right)^2$, which has a global minimum at $\mathbf{1}$.
We use a step size of $8{\cdot}10^{-5}$ for SGD and $0.9$ for the binary gradient optimizers signSGD,  (D) and  (G).
<ref> demonstrates the convergence of various algorithms.
As the theoretical convergence rates are identical, we see a similar trend for all methods.
For this toy problem, the best convergence is achieved by  where the weight updates are computed on the Gurobi optimizer, while the variant with D-Wave QA optimizer is performing as well as signSGD.
This empirical evidence (i) justifies the validity of quantum hardware in obtaining good-quality solutions, and (ii) validates our continuous update rule.

§.§ Training binary logistic regression

Before leaping to neural networks, we evaluate the performance of our algorithm on the well-studied binary logistic regression, a logistic regression model where the target variable is binary.
We synthesize a toy dataset to show empirical convergence. First, we sample two, two-dimensional linearly separable blobs, label them with $0$ and $1$, accordingly and lift them into a $3D$ space by appending the $\mathbf{1}$ intercept.
Next, on this data, we fit the logistic regression model using our algorithm as well as the baselines.
Our model consists of a binary linear layer with an input dimension of $3$ and an output dimension of $1$, followed by a

Binary logistic regression.

sigmoid activation function to perform binary classification. It is trained with a binary cross entropy (BCE) loss.
We use the following learning rates:
$0.05$ for  (Gurobi) and BC (signSGD), $0.01$ for  (D-Wave) and $5{\cdot}10^{-5}$ for BC (SGD).
In <ref>, we observe that after the initial fitting phase, our optimizer is consistently the best-performing one when the binary projection is computed via Gurobi.
This indicates two notable points: (i) our algorithm benefits from solving a QUBO to find the best fitting binary gradient approximation, compared to using $\sign(\cdot)$, and (ii) we are dependent on the quality of projections, , as D-Wave solution has lower quality, our convergence is slower.

§.§ Training of binary MLPs

We now present training of binary MLPs for binary classification on two datasets: UCI Adult [53] and MNIST [25]. We experimentally validate the effectiveness of our approach in training binary multi-layer perceptions by (i) demonstrating that  converges and (ii) showing how well the MLPs perform when trained with the different optimizers.

§.§.§ Training on UCI Adult

The UCI Adult [53] is a binary classification benchmark containing $123$ binary features, $1,605$ labelled training data and $30,956$ testing data. We now analyze the convergence and training accuracy on certain subsets of this dataset. Models in this setup are trained with the BCE loss.

Mean accuracy over five runs. Left: two- layer setup; Right: five-layer setup. We compare  against ProxQuant, BinaryConnect (BC) with signSGD & SGD optimizer.

Convergence of .
We now investigate the convergence of a two-layer and a five-layer BNN with $15$ input dimensions, $10$ hidden, and one output dimension.
As input, we use a UCI subset of 15 random features and 32 batches with a batch size of 16.
<ref> demonstrates the convergence of our algorithm in comparison to our baselines. As the training is dominated by noise, we show here an average of over five runs.
The convergence of  is similar to BC. Yet, while ProxQuant initially converges faster in the two-layer setup (due to the real weights that gradually convert to binary), it fails to train in the five-layer setup. In the two-layer setup, BC-SGD seems to have a slight convergence advantage. However, in the five-layer network, we observe a tie between the BC versions and . While the theoretical convergence rates are identical, ours is the only algorithm with the ability to leverage QCs.

Training accuracy on a subset of the UCI Adult dataset for binary classification.

Training accuracy.
We now use three manually selected features and 30 batches with a batch size of 20 to train a two-layer neural network with the input dimension of $3$, hidden dimension of $5$ and the output of $2$.
<ref> reports the training accuracy where we observe that  (G) outperforms the baselines as well as our D-Wave optimiser,  (D), which nevertheless performs better than or on par with the binary signSGD.

§.§.§ Training on MNIST

max width=

ProxQuant    BC     BC

SGD      signSGD     (Gurobi)     (D-Wave)

0/2      0.65    0.64     0.71     0.66       0.62

1/2    0.67     0.72     0.66     0.73     0.70

1/7     0.64    0.74     0.68     0.75     0.74

The accuracy of binary classification on MNIST. The first column contains the digits used in the experiment.

We now turn to a famous problem of digit classification on MNIST dataset [41].
To ensure that our features are compact and binary we extract handcrafted features of MNIST letters as follows:
(i) we extract key points via Monti  [48], (ii) we sample 16 lines, which all run through the centroid of the image
and (iii) for each line, we extract a descriptor by counting the number of key points falling in the positive and negative half-plane. If there are more key points on the positive half-plane of line $i$ we set $s_i=1$ and otherwise to $s_i=-1$. For a given image, this amounts to a feature vector $\mathbf{s}\in\{\pm 1\}^{16}$.
In total, we use a train/test split of $500/3000$ images.
We then adopt a three-layer MLP, with an input dimension of $16$, a hidden dimension of $4$, and an output of $2$ to binary-classify the digits of $0$, $1$, $2$ and $7$ into two classes via negative log-likelihood (NLL) loss. The accuracies attained are reported in <ref>.
On this setup,  (D) shows to be on par with the state of the art, while our  (G) can achieve the best accuracy. With the advancements in quantum hardware, we expect  (D) to surpass its classical counterpart.

§.§ Training of binary GCNNs

Finally, we present how our algorithm performs on graph data.
We start with tests on the small-scale example of the Karateclub social network graph [77].
In addition, we report the performance of our method on the citation graphs of Cora [46] and Pubmed [50] benchmarks.
We start by providing our experimental configuration.

§.§.§ Datasets and Implementation Details

Karate club
We are using the Karate club social network graph [77] with handcrafted features.
The feature vectors one-hot encode the class and thus form an identity matrix.
We reduce those vectors by changing the one hot encoding to a binary encoding, e.g., $0010 \rightarrow 10$.
This reduces the dimensions from $34$ down to $6$.
As a train/test split, we use 20 nodes selected to maintain a balance across all four labels with an equal number of nodes assigned to each label, and we reserve 14 nodes for testing.

Cora and Pubmed
We use training/validation/test splits in accordance with Yang  [72]: 140/500/1000 for Cora and 60/500/1000 for Pubmed, where the remaining nodes are treated as unlabelled data.
During training, we use minibatches comprising $1/16^{\mathrm{th}}$ of the training data, from which we draw a single sample per epoch. To reduce the dimensionality of the convolution layers (to $10$ in this case), we use a two-layer MLP. This encoder-MLP is in all cases trained with Adam. Nevertheless, the resulting layer-wise QUBO problem, even in the low dimensions, is too large for D-Wave Advantage 6.1. Hence, we use the Gurobi solver to demonstrate the efficacy of our algorithm, .
On all three setups, we employ the NLL loss.

Graph classification: We report mean test accuracy over five runs for Karate club [77] (left), Cora [46] (middle) and Pubmed [50] (right) datasets.

$N_\ell = 1$ neuron

Close-up of (a)

$N_\ell = 2$ neurons

Close-up of (c)

The eigenvalues of the Hamiltonian of the QUBO in <ref> while training a binary MLP on the adult dataset. Those eigenvalues are plotted as a function of annealing time ($t/T$) for a linear layer with one to four neurons and batch size 1. The red bar represents the eigenvalue gap between the ground level and the first excited level that does not evolve into the ground state.

§.§.§ Discussion

Our results are plotted in <ref>.
First to notice is that on all three datasets  (G) outperforms the baselines on the average.
SGD is already known to be a suboptimal optimiser when learning on graphs [34].
We see that BinaryConnect (BC) also inherits this property of SGD when trained with SGD.
The optimality of our weight update scheme, therefore, improves the effectiveness of training significantly compared to the BC-SGD.
While BC (signSGD) and ProxQuant exhibit faster convergence compared to , our approach achieves superior overall accuracy across all datasets. This once again underscores the inherent advantage of employing a QUBO framework for attaining optimal binary weight updates.
Furthermore, the remaining techniques lack quantum deployability, making  the pioneering quantum-annealing-executable algorithm, as far as our knowledge extends. While it is foreseeable that our algorithm encounters limitations imposed either by classical QUBO solvers or quantum hardware, it consistently outperforms classical specialised algorithms such as BinaryConnect and ProxQuant. Further, we expect our algorithm to inherit the advances in quantum computing and naturally improve over time.

§.§ Eigenvalue Gap

We now examine, via simulation, the eigenvalue gap corresponding to the QUBO in <ref> for small MLPs. <ref> plots the eigenvalues (blue lines) of the time-varying Hamiltonian as a function of $t/T$, for a single binary linear layer $l$ with up to four neurons and batch size 1 (besides these parameters, the size of $H$, and hence the number of distinct eigenvalues, depend on other factors; more on the feasibility of simulation below). The smallest eigenvalue gap is highlighted by a red bar in the plots. With the same annealing time, a larger eigenvalue gap yields a final state closer to the ground state and, thus, a higher-quality solution.

The results show that, as expected, the eigenvalue gap shrinks gradually with increasing problem size. While four data points are not enough to extrapolate a trend, the exercise shows how the performance of quantum annealing depends on the hyperparameters. In particular, the observations suggest keeping the number of neurons of each layer
small, which argues for updating each layer partially.

§ CONCLUSION

We introduced , a scalable hybrid quantum-classical optimiser designed for incremental training of binary neural networks on real quantum hardware. By capitalising on the nature of the optimisation problem associated with incremental binary weight update, we have developed a projection procedure that admits a QUBO form effectively solvable via QA. Under practically acceptable assumptions we have also provided convergence proofs validating the reliability of our algorithm. Our results demonstrate that  achieves comparable or superior performance compared to existing state of the art, while remaining adaptable to the advances in quantum computation.

The problems we can solve on quantum annealers are size-limited by the current capacity of the quantum processing units.
Nevertheless, we were able to optimise the network weights in
all our experiments.
Moreover, due to the rapid improvements in quantum hardware, we are optimistic that our algorithm can impact practical applications already in close future.
At times, we have found the classical SA solutions to be superior to QAs. We expect this temporary advantage to flip in favor of QAs with improvements of the D-Wave hardware or QAs in general.

Future work
Promising avenues involve (i) quantum noise reduction techniques to enhance the D-Wave solution, (ii) adding momentum to our optimiser and (iii) a divide and conquer approach to train larger networks.

Acknowledgements Tolga Birdal and Vladislav Golyanik acknowledge compute resource support from D-Wave. This work was supported by the AWS Cloud Credit for Research program.

[1]
Steven H Adachi and Maxwell P Henderson.
Application of quantum annealing to training of deep neural networks.
arXiv preprint arXiv:1510.06356, 2015.

[2]
Federica Arrigoni, Willi Menapace, Marcel Seelbach Benkner, Elisa Ricci, and
Vladislav Golyanik.
Quantum motion segmentation.
In European Conference on Computer Vision (ECCV), 2022.

[3]
Davis Arthur, Lauren Pusey-Nazzaro, et al.
Qubo formulations for training machine learning models.
Scientific reports, 11(1):1–10, 2021.

[4]
Mayowa Ayodele, Richard Allmendinger, Manuel López-Ibáñez, and
Matthieu Parizy.
Multi-objective qubo solver: bi-objective quadratic assignment
In Proceedings of the Genetic and Evolutionary Computation
Conference, pages 467–475, 2022.

[5]
Mehdi Bahri, Gaétan Bahl, and Stefanos Zafeiriou.
Binary graph neural networks.
In CVPR, pages 9492–9501, 2021.

[6]
Yu Bai, Yu-Xiang Wang, and Edo Liberty.
Proxquant: Quantized neural networks via proximal operators.
arXiv preprint arXiv:1810.00861, 2018.

[7]
Kerstin Beer, Dmytro Bondarenko, Terry Farrelly, Tobias J Osborne, Robert
Salzmann, Daniel Scheiermann, and Ramona Wolf.
Training deep quantum neural networks.
Nature communications, 11(1):1–6, 2020.

[8]
Yoshua Bengio, Nicholas Léonard, and Aaron Courville.
Estimating or propagating gradients through stochastic neurons for
conditional computation.
arXiv preprint arXiv:1308.3432, 2013.

[9]
Harshil Bhatia, Edith Tretschk, Zorah Lähner, Marcel Benkner, Michael Möller,
Christian Theobalt, and Vladislav Golyanik.
Ccuantumm: Cycle-consistent quantum-hybrid matching of multiple
In Computer Vision and Pattern Recognition (CVPR), 2023.

[10]
Jacob Biamonte, Peter Wittek, Nicola Pancotti, Patrick Rebentrost, Nathan
Wiebe, and Seth Lloyd.
Quantum machine learning.
Nature, 549(7671):195–202, 2017.

[11]
Tolga Birdal, Vladislav Golyanik, Christian Theobalt, and Leonidas Guibas.
Quantum permutation synchronization.
In CVPR, 2021.

[12]
Kelly Boothby, Paul Bunyk, Jack Raymond, and Aidan Roy.
Next-generation topology of d-wave quantum processors.
arXiv e-prints, 2020.

[13]
Nicolas Boumal.
An introduction to optimization on smooth manifolds.
Cambridge University Press, 2023.

[14]
Michael Broughton, Guillaume Verdon, Trevor McCourt, Antonio J Martinez,
Jae Hyeon Yoo, Sergei V Isakov, Philip Massey, Ramin Halavati, Murphy Yuezhen
Niu, Alexander Zlokapa, et al.
Tensorflow quantum: A software framework for quantum machine
arXiv preprint arXiv:2003.02989, 2020.

[15]
Adrian Bulat and Georgios Tzimiropoulos.
Xnor-net++: Improved binary neural networks.
arXiv preprint arXiv:1909.13863, 2019.

[16]
Zhaowei Cai, Xiaodong He, Jian Sun, and Nuno Vasconcelos.
Deep learning with low precision by half-wave gaussian quantization.
In CVPR, pages 5918–5926, 2017.

[17]
Tat-Jun Chin, David Suter, Shin-Fang Chng, and James Quach.
Quantum robust fitting.
In ACCV, 2020.

[18]
Jungwook Choi, Zhuo Wang, Swagath Venkataramani, Pierce I-Jen Chuang,
Vijayalakshmi Srinivasan, and Kailash Gopalakrishnan.
Pact: Parameterized clipping activation for quantized neural
arXiv preprint arXiv:1805.06085, 2018.

[19]
Ronan Collobert, Jason Weston, Léon Bottou, Michael Karlen, Koray
Kavukcuoglu, and Pavel Kuksa.
Natural language processing (almost) from scratch.
Journal of machine learning research, 12(ARTICLE):2493–2537,

[20]
Iris Cong, Soonwon Choi, and Mikhail D Lukin.
Quantum convolutional neural networks.
Nature Physics, 15:1273–1278, 2019.

[21]
Matthieu Courbariaux, Yoshua Bengio, and Jean-Pierre David.
Binaryconnect: Training deep neural networks with binary weights
during propagations.
In NIPS, pages 3123–3131, 2015.

[22]
<https://docs.dwavesys.com/docs/latest/doc_handbook.html>, 2021.

[23]
<https://docs.dwavesys.com/docs/latest/c_gs_4.html>, 2021.

[24]
Yann N. Dauphin, Harm de Vries, Junyoung Chung, and Yoshua Bengio.
Rmsprop and equilibrated adaptive learning rates for non-convex
CoRR, 2015.

[25]
Li Deng.
The mnist database of handwritten digit images for machine learning
IEEE Signal Processing Magazine, 29(6):141–142, 2012.

[26]
David P DiVincenzo.
Quantum computation.
Science, 270(5234):255–261, 1995.

[27]
Vivek Dixit, Raja Selvarajan, Muhammad A Alam, Travis S Humble, and Sabre Kais.
Training restricted boltzmann machines with a d-wave quantum
Front. Phys. 9: 589626. doi: 10.3389/fphy, 2021.

[28]
Matteo Farina, Luca Magri, Willi Menapace, Elisa Ricci, Vladislav Golyanik, and
Federica Arrigoni.
Quantum multi-model fitting.
In Computer Vision and Pattern Recognition (CVPR), 2023.

[29]
Siddhant Garg and Goutham Ramakrishnan.
Advances in quantum deep learning: An overview.
arXiv preprint arXiv:2005.04316, 2020.

[30]
Amir Gholami, Sehoon Kim, Zhen Dong, Zhewei Yao, Michael W Mahoney, and Kurt
A survey of quantization methods for efficient neural network
arXiv preprint arXiv:2103.13630, 2021.

[31]
Lu Hou, Quanming Yao, and James T Kwok.
Loss-aware binarization of deep networks.
arXiv preprint arXiv:1611.01600, 2016.

[32]
Linyong Huang, Zhe Zhang, Zhaoyang Du, Shuangchen Li, Hongzhong Zheng, Yuan
Xie, and Nianxiong Tan.
Epquant: A graph neural network compression approach based on product
Neurocomputing, 503:49–61, 2022.

[33]
Itay Hubara, Matthieu Courbariaux, Daniel Soudry, Ran El-Yaniv, and Yoshua
Binarized neural networks.
NIPS, 29, 2016.

[34]
Mohammad Rasool Izadi, Yihao Fang, Robert Stevenson, and Lizhen Lin.
Optimization of graph neural networks with natural gradient descent.
In 2020 IEEE international conference on big data (big data),
pages 171–179. IEEE, 2020.

[35]
Yongcheng Jing, Yiding Yang, Xinchao Wang, Mingli Song, and Dacheng Tao.
Meta-aggregator: learning to aggregate for 1-bit graph neural
In CVPR, pages 5301–5310, 2021.

[36]
Iordanis Kerenidis, Jonas Landman, and Anupam Prakash.
Quantum algorithms for deep convolutional neural networks.
arXiv preprint arXiv:1911.01117, 2019.

[37]
B Kerstin, B Dmytro, F Terry, O Tobias, S Robert, and W Ramona.
Efficient learning for deep quantum neural networks.
Nature, 2019.

[38]
Diederik P Kingma and Jimmy Ba.
Adam: A method for stochastic optimization.
arXiv preprint arXiv:1412.6980, 2014.

[39]
Thomas N Kipf and Max Welling.
Semi-supervised classification with graph convolutional networks.
arXiv preprint arXiv:1609.02907, 2016.

[40]
Alex Krizhevsky, Ilya Sutskever, and Geoffrey E Hinton.
Imagenet classification with deep convolutional neural networks.
Communications of the ACM, 60(6):84–90, 2017.

[41]
Yann LeCun and Corinna Cortes.
MNIST handwritten digit database.
http://yann.lecun.com/exdb/mnist/, 2010.

[42]
M. Lichman.
Uci machine learning repository.
<https://archive.ics.uci.edu/ml>, 2013.

[43]
Zechun Liu, Zhiqiang Shen, Marios Savvides, and Kwang-Ting Cheng.
Reactnet: Towards precise binary neural network with generalized
activation functions.
In ECCV, pages 143–159. Springer, 2020.

[44]
Zechun Liu, Baoyuan Wu, Wenhan Luo, Xin Yang, Wei Liu, and Kwang-Ting Cheng.
Bi-real net: Enhancing the performance of 1-bit cnns with improved
representational capability and advanced training algorithm.
In ECCV, pages 722–737, 2018.

[45]
Andrew Lucas.
Ising formulations of many np problems.
Frontiers in physics, 2:5, 2014.

[46]
Andrew Kachites McCallum, Kamal Nigam, Jason Rennie, and Kristie Seymore.
Automating the construction of internet portals with machine
Information Retrieval, 3:127–163, 2000.

[47]
Timothy M McCormick, Zipporah Klain, Ian Herbert, Anthony M Charles, R Blair
Angle, Bryan R Osborn, and Roy L Streit.
Multiple target tracking and filtering using bayesian diabatic
quantum annealing.
arXiv preprint arXiv:2209.00615, 2022.

[48]
Federico Monti, Davide Boscaini, Jonathan Masci, Emanuele Rodola, Jan Svoboda,
and Michael M Bronstein.
Geometric deep learning on graphs and manifolds using mixture model
In CVPR, pages 5115–5124, 2017.

[49]
Markus Nagel, Marios Fournarakis, Rana Ali Amjad, Yelysei Bondarenko, Mart van
Baalen, and Tijmen Blankevoort.
A white paper on neural network quantization.
arXiv preprint arXiv:2106.08295, 2021.

[50]
Galileo Namata, Ben London, Lise Getoor, Bert Huang, and U Edu.
Query-driven active surveying for collective classification.
In 10th international workshop on mining and learning with
graphs, volume 8, page 1, 2012.

[51]
Michael A Nielsen and Isaac Chuang.
Quantum computation and quantum information, 2002.

[52]
Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory
Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban
Desmaison, Andreas Kopf, Edward Yang, Zachary DeVito, Martin Raison, Alykhan
Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith
Pytorch: An imperative style, high-performance deep learning library.
In NeurIPS, pages 8024–8035. Curran Associates, Inc., 2019.

[53]
John Platt.
Sequential minimal optimization: A fast algorithm for training
support vector machines.
Technical Report MSR-TR-98-14, Microsoft, April 1998.

[54]
Haotong Qin, Ruihao Gong, Xianglong Liu, Xiao Bai, Jingkuan Song, and Nicu
Binary neural networks: A survey.
Pattern Recognition, 105:107281, 2020.

[55]
Haotong Qin, Ruihao Gong, Xianglong Liu, Xiao Bai, Jingkuan Song, and Nicu
Binary neural networks: A survey.
Pattern Recognition, 105, 2020.
arXiv: 2004.03333.

[56]
Mohammad Rastegari, Vicente Ordonez, Joseph Redmon, and Ali Farhadi.
Xnor-net: Imagenet classification using binary convolutional neural
In European conference on computer vision, pages 525–542.
Springer, 2016.

[57]
Lakshika Rathi, Edith Tretschk, Christian Theobalt, Rishabh Dabral, and
Vladislav Golyanik.
3D-QAE: Fully quantum auto-encoding of 3d point clouds.
In BMVC, 2023.

[58]
Mher Safaryan and Peter Richtárik.
Stochastic sign descent methods: New algorithms and better theory.
In International Conference on Machine Learning, pages
9224–9234. PMLR, 2021.

[59]
Michele Sasdelli and Tat-Jun Chin.
Quantum annealing formulation for binary neural networks.
In 2021 Digital Image Computing: Techniques and Applications
(DICTA), pages 1–10. IEEE, 2021.

[60]
Maria Schuld, Ilya Sinayskiy, and Francesco Petruccione.
The quest for a quantum neural network.
Quantum Information Processing, 13(11):2567–2586, 2014.

[61]
Maria Schuld, Ilya Sinayskiy, and Francesco Petruccione.
An introduction to quantum machine learning.
Contemporary Physics, 56(2):172–185, 2015.

[62]
Marcel Seelbach Benkner, Vladislav Golyanik, Christian Theobalt, and
Michael Moeller.
Adiabatic quantum graph matching with permutation matrix constraints.
In International Conference on 3D Vision (3DV), 2020.

[63]
Wei Tang, Gang Hua, and Liang Wang.
How to train a compact binary neural network with high accuracy?
In Thirty-First AAAI conference on artificial intelligence,

[64]
Tijmen Tieleman and Geoffrey Hinton.
Lecture 6.5-rmsprop: Divide the gradient by a running average of its
recent magnitude, 2012.

[65]
Jirí Tumpach, Jan Kalina, and Martin Holena.
A comparison of regularization techniques for shallow neural
networkstrained on small datasets.
In ITAT, pages 94–103, 2021.

[66]
Guillaume Verdon, Trevor McCourt, Enxhell Luzhnica, Vikash Singh, Stefan
Leichenauer, and Jack Hidary.
Quantum graph neural networks.
arXiv preprint arXiv:1909.12264, 2019.

[67]
Hanchen Wang, Defu Lian, Ying Zhang, Lu Qin, Xiangjian He, Yiguang Lin, and
Xuemin Lin.
Binarized graph neural network.
World Wide Web, 24:825–848, 2021.

[68]
Junfu Wang, Yunhong Wang, Zhen Yang, Liang Yang, and Yuanfang Guo.
Bi-gcn: Binary graph convolutional network.
In CVPR, pages 1561–1570, 2021.

[69]
Shaopeng Wei and Yu Zhao.
Graph learning: A comprehensive survey and future directions.
arXiv preprint arXiv:2212.08966, 2022.

[70]
Nathan Wiebe, Ashish Kapoor, and Krysta M Svore.
Quantum deep learning.
arXiv preprint arXiv:1412.3489, 2014.

[71]
Peter Wittek.
Quantum machine learning: what quantum computing means to data
Academic Press, 2014.

[72]
Zhilin Yang, William Cohen, and Ruslan Salakhudinov.
Revisiting semi-supervised learning with graph embeddings.
In ICML, pages 40–48. PMLR, 2016.

[73]
Kai-Lang Yao and Wu-Jun Li.
Full-precision free binary graph neural networks, 2022.

[74]
Yang You, Jing Li, Sashank Reddi, Jonathan Hseu, Sanjiv Kumar, Srinadh
Bhojanapalli, Xiaodan Song, James Demmel, Kurt Keutzer, and Cho-Jui Hsieh.
Large batch optimization for deep learning: Training bert in 76
arXiv preprint arXiv:1904.00962, 2019.

[75]
Chunyu Yuan and Sos S Agaian.
A comprehensive review of binary neural network.
arXiv preprint arXiv:2110.06804, 2021.

[76]
Alp Yurtsever, Tolga Birdal, and Vladislav Golyanik.
Q-fw: A hybrid classical-quantum frank-wolfe for quadratic binary
arXiv preprint arXiv:2203.12633, 2022.

[77]
Wayne W Zachary.
An information flow model for conflict and fission in small groups.
Journal of anthropological research, 33(4):452–473, 1977.

[78]
Jan-Nico Zaech, Alexander Liniger, Martin Danelljan, Dengxin Dai, and Luc
Van Gool.
Adiabatic quantum computing for multi object tracking.
In CVPR, 2022.

[79]
Dongqing Zhang, Jiaolong Yang, Dongqiangzi Ye, and Gang Hua.
Lq-nets: Learned quantization for highly accurate and compact deep
neural networks.
In ECCV, pages 365–382, 2018.

[80]
Ziwei Zhang, Peng Cui, and Wenwu Zhu.
Deep learning on graphs: A survey.
IEEE Transactions on Knowledge and Data Engineering,
34(1):249–270, 2020.

[81]
Shuchang Zhou, Yuxin Wu, Zekun Ni, Xinyu Zhou, He Wen, and Yuheng Zou.
Dorefa-net: Training low bitwidth convolutional neural networks with
low bitwidth gradients.
arXiv preprint arXiv:1606.06160, 2016.

[82]
Zeyu Zhu, Fanrong Li, Zitao Mo, Qinghao Hu, Gang Li, Zejian Liu, Xiaoyao Liang,
and Jian Cheng.
A2q: Aggregation-aware quantization for graph neural networks.
arXiv preprint arXiv:2302.00193, 2023.

§ APPENDIX

This documents supplements our paper with the additional material referred in the main text. This includes the proofs of the propositions (specifically the Euclidean to binary gradient, $\rho-$norm, the convergence of SBGD and P-SBGD, and the QUBO derivation), supporting evidence of the directionality and additional experimental evaluations, including visualizations of the Rosenbrock function, the training of real-valued MLPs, and analysis of the D-Wave QUBO problems. In addition, we provide an exploration of the behaviour of the learning rate and batch size and furnish extra information on related work, implementation specifics, and function definitions.

§ THEORETICAL RESULTS

§.§ Approximating Euclidean gradients with binary gradients

We start by re-phrasing <ref> defined in the main paper as a proposition proving the equality for the continuous relaxation of the binary map to the gradient.

Consider a matrix $\bZ^t\in\R^{n\times m}$ whose $i^\mathrm{th}$ column is $\R^n\ni\bz_i^t = \bj^t_i / \| \bj_i^t\|^2 $ where $\bj_i^t = \frac{\partial y_i}{\partial \x}|_{\x^t}$ and $\y=f(\x)$. Note that, $\bj_i^t$ denotes the $i^\mathrm{th}$ column of the Jacobian map $\R^{n\times m}\ni\J^t:= \frac{\partial \y}{\partial \x}$ evaluated at $\x^t$.
Then, using <ref> to map the gradient w.r.t. the output $\y$ allows us to get the vector approximating the gradient w.r.t. the parameters $\x$:

\begin{align}
\hat{\bPi}_{\bZ^t}(\tilde{\nabla}_{\y}E_f |_{\x^t}) = \tilde{\nabla}_{\x}E_f |_{\x^t}.
\label{eq:egrad2bgradRelaxed}
\end{align}

In the case binary regime, we resort to its approximation which involves our binary projection:

\begin{align}
\bPi_{\bZ^t}(\tilde{\nabla}_{\y}E_f |_{\x^t}) \approx \argmin_{\bv\in\{-1,1\}^{n}} \|\tilde{\nabla}_{\x}E_f |_{\x^t} -\bv\|_2^2.
\end{align}

We now give the proof for the continuous map given in <ref>, where we optimise over $\R^{n}$. We now expand <ref>

\begin{align}
\frac{\partial E}{\partial \x}= \argmin_{\bg \in \R^{n}}\sum_{i = 1}^m|\frac{\partial E}{\partial y_i} - \bg^\top \z_i^t |^2.
\label{eq:objective_nn}
\end{align}

This statement is true if every term in the sum is $0$. If this holds, we can replace $\bg$ by $\frac{\partial E}{\partial \x}$ and seek to have

\begin{align}
\frac{\partial E}{\partial y_i} - \frac{\partial E}{\partial \x}^\top\bz_i^t  = 0.
\label{eq:statement_to_show}
\end{align}

Assuming differentiability and compositionality we can use the chain rule to write

\begin{equation}\frac{\partial E}{\partial \x} = \frac{\partial E}{\partial y_i} \cdot \frac{\partial y_i}{\partial \x}.
\label{eq:splitGradient}
\end{equation}

By <ref>, and the definition of $\z_i^t$ we can reformulate the left hand side of <ref> as follows:

\begin{align}
\frac{\partial E}{\partial y_i} - \left(\frac{\partial E}{\partial \x}\right)^\top  \left(\frac{\partial y_i}{\partial \x} / \|\frac{\partial y_i}{\partial \x} \|^2\right) &=\frac{\partial E}{\partial y_i} - \frac{\partial E}{\partial y_i}\left(\frac{\partial y_i}{\partial \x}\right)^\top \left(\frac{\partial y_i}{\partial \x} / \|\frac{\partial y_i}{\partial \x} \|^2\right) \\
&= \frac{\partial E}{\partial y_i}\left(1 -  \left(\frac{\partial y_i}{\partial \x} / \|\frac{\partial y_i}{\partial \x} \|\right)^\top \left(\frac{\partial y_i}{\partial \x} / \|\frac{\partial y_i}{\partial \x} \|\right)\right)
\label{eq:start_of_matrixes}\\
&=\frac{\partial E}{\partial y_i}(1 - 1) = 0 \qquad \mathrm{for\,all} \,\, i.\label{eq:jacmap}
\end{align}

Alternatively, we can arrive at the same proof by using the fact that Jacobians map tangent vectors, , local isomorphisms between the tangent spaces of input and output points. Finally, by satisfying <ref> for all $i$, we see that in the continuous space, <ref> is fulfilled.

We can use this mapping to calculate binary updates also for function parameters, which are matrices and not just vectors.
This can easily be done by vectorising the parameter matrix into $m$ vectors and computing $m$ times the mapping of the gradient of a parameter vector.

§.§ Proof that the $\rho$-norm is a norm

We first restate the proposition <ref> of the main paper before we provide the full proof.

If CDP holds, then <ref> is a norm in the mathematical sense,  fulfills subadditivity, homogeneity and positive definiteness.

Let $\g(\x) = \nabla_{\x} f(\x)$. We now recall the definition of $\rho$-norm:

\begin{align}
\|\bg(\x)\|_\rho:=\sum_{i=1}^d\left(2 \rho_i(\x)-1\right)\left|g_i(\x)\right|
\label{def:rhonorm_appendix}
\end{align}

and prove it by ensuring that the properties of a norm hold: subadditivity, absolute homogeneity and positive definiteness.

Subadditivity Let $\mathbf{g},\mathbf{h}$ be two gradients corresponding with the same CDP. Hence, they both and the sum of the two share the same probability function $\rho$. We then have:

\begin{align}
\|\bg(\x) + \mathbf{h}(\x) \|_\rho&=  \sum_{i=1}^d\left(2 \rho_i(\x)-1\right)\left|g_i(\x) + h_i(\x)\right| \\
&\leq \sum_{i=1}^d\left(2 \rho_i(\x)-1\right)(\left|g_i(\x)\right| +\left| h_i(\x)\right|) \\
&=  \sum_{i=1}^d\left(2 \rho_i(\x)-1\right)\left|g_i(\x)\right| +\sum_{i=1}^d\left(2 \rho_i(\x)-1\right)\left| h_i(\x)\right| \\
&=  \|\bg(\x) \|_\rho  +  \| \mathbf{h}(\x) \|_\rho.
\end{align}

Absolute homogeneity

\begin{align*}
\|s\bg(\x)  \|_\rho &=  \sum_{i=1}^d\left(2 \rho_i(\x)-1\right)\left|sg_i(\x)\right| \\
&\leq \sum_{i=1}^d\left(2 \rho_i(\x)-1\right)|s|\left|g_i(\x)\right| \\
&= |s|\sum_{i=1}^d\left(2 \rho_i(\x)-1\right)\left|g_i(\x)\right| \\
&= |s|\|\bg(\x)  \|_\rho.
\end{align*}

Positive definiteness
As the CDP holds: $\forall i: \left(2 \rho_i(\x)-1\right) > 0$. Hence $ \sum_{i=1}^d\left(2 \rho_i(\x)-1\right)\left|g_i(\x)\right|$ is only equal to 0, when : $\forall i: g_i(\x) = 0$. This is only the case for $\bg = 0$. Hence the positive definiteness is shown as well.

§.§ Proof of convergence of SBGD

Let us restate our theorem before providing our proof:

Assuming H1 [Lower bound], and H2 [Co-Directionality Probability (CDP)] and using a step size of $\alpha_t=\alpha_0 / \sqrt{t+1}$,  guarantees a convergence rate of $\mathcal{O}({1}/{\sqrt{T}})$:

\begin{align}
\min _{0 \leq t<T} \mathbb{E}\left\|\nabla_{\x} E_f\left(\x^t\right)\right\|_\rho \leq \frac{E_f\left(\x^0\right)-E_f^*}{\alpha_0 \sqrt{T}}+\frac{3 \alpha_0 n \bar{K}}{2} \frac{\log T}{\sqrt{T}}.
\end{align}

If $\alpha_t \equiv \alpha>0$, we get a convergence of $\mathcal{O} ({1}/{T})$ convergence to a neighbourhood:

\begin{align}
\frac{1}{T} \sum_{t=0}^{T-1} \mathbb{E}\left\|\nabla_{\x} E_f\left(\x^t\right)\right\|_\rho \leq \frac{E_f\left(\x^0\right)-E_f^*}{\alpha T}+\frac{\alpha n \bar{K}}{2}.
\end{align}

Our proof is mostly based upon Safaryan and Richtárik [58]. We will start by using the Lipschitz smoothness assumption to obtain

\begin{align}
E_f\left(\x^{t+1}\right)& =E_f\left(\x^t-\alpha_t \bPi_{\bZ^t}( \widetilde{\nabla}_{\y} E_f(\x^t))\right) \\
& \leq E_f\left(\x^t\right)-\left\langle \bg_t, \alpha_t \bPi_{\bZ^t}( \widetilde{\nabla}_{\y}E_f(\x^t))\right\rangle+\sum_{i=1}^n \frac{K_i}{2}\left(\alpha_t \bPi_{\bZ^t}( \widetilde{\nabla}_{\y}E_f(\x^t))_i\right)^2 \\
& =E_f\left(\x^t\right)-\alpha_t\left\langle g_t, \bPi_{\bZ^t}( \widetilde{\nabla}_{\y}E_f(\x^t))\right\rangle+\frac{n \bar{K}}{2} \alpha_t^2
\end{align}

where $\bg_t=\bg\left(\x^t\right) = \nabla_{\x} E_f\left(\x^t\right)$. Taking conditional expectation given current iteration $\x^t$ gives

\begin{align}
\begin{split}
&\mathbb{E}\left[E_f\left(\x^{t+1}\right) \mid
\x^t\right] \leq E_f\left(\x^t\right)-\alpha_t
\mathbb{E}\left[\left\langle \bg_t,
\bPi_{\bZ^t}
(\widetilde{\nabla}_\y E_f(\x^t))
\right\rangle\right]+\frac{n \bar{K}}{2} \alpha_t^2.
\label{eq:main_inequality}
\end{split}
\end{align}

We will now show that the equality $\mathbb{E}\left[\left\langle \nabla_{\bx}E_f(\x^t), \bPi_{\bZ^t}( \widetilde{\nabla}_{\y}E_f(\x^t))\right\rangle\right]=\left\|\nabla_{\bx}E_f(\x^t)\right\|_\rho$ holds.

\begin{align}
\mathbb{E}\left[\left\langle \bg_t, \bPi_{\bZ^t}
(\widetilde{\nabla}_\y E_f(\x^t))\right\rangle\right]
=&\left\langle \bg_t, \mathbb{E}\left[\bPi_{\bZ^t}
(\widetilde{\nabla}_\y E_f(\x^t))\right]\right\rangle \\
=&\sum_{i=1}^n g_{t, i}  \mathbb{E}\left[\bPi_{\bZ^t}
(\widetilde{\nabla}_\y E_f(\x^t))_i\right] \\
=&\sum_{\substack{1 \leq i \leq n \\
g_{t, i} \neq 0}} g_{t, i}  \mathbb{E}\left[\bPi_{\bZ^t}
(\widetilde{\nabla}_\y E_f(\x^t))_i\right]. \label{eq:last_eq}
\end{align}

In the following, we can now use the CDP hypothesis and definition. This transforms <ref> to:

\begin{align}
&\sum_{\substack{1 \leq i \leq n \\
g_{t, i} \neq 0}} g_{t, i}\left(\rho_i\left(\x^t\right) \operatorname{sign} g_{t, i}+\left(1-\rho_i\left(\x^t\right)\right)\left(-\operatorname{sign} g_{t, i}\right)\right) \\
=&\sum_{\substack{1 \leq i \leq n \\
g_{t, i} \neq 0}}\left(2 \rho_i\left(\x^t\right)-1\right)\left|g_{t, i}\right|\\
=&\sum_{i=1}^n\left(2 \rho_i\left(\x^t\right)-1\right)\left|g_{t, i}\right|=\left\|\bg_t\right\|_\rho.
\end{align}

Hence the equation $\mathbb{E}\left[\left\langle \nabla_{\bx}E_f(\x^t), \bPi_{\bZ^t}( \widetilde{\nabla}_{\y}E_f(\x^t))\right\rangle\right]=\left\|\nabla_{\bx}E_f(\x^t)\right\|_\rho$
is shown.
We can now plug this result into <ref>, taking the expectations on both sides and rearranging the individual terms we get:

\begin{align}
\mathbb{E}\left\|\bg_t\right\|_\rho \leq \frac{\mathbb{E}\left[E_f\left(\x^t\right)\right]-\mathbb{E}\left[E_f\left(\x^{t+1}\right)\right]}{\alpha_t}+\frac{n \bar{K}}{2} \alpha_t.
\end{align}

\begin{align}
\sum_{t=0}^{T-1} \alpha_t \mathbb{E}\left\|\bg_t\right\|_\rho \leq\left(E_f\left(\x^0\right)-E_f^*\right)+\frac{n \bar{K}}{2} \sum_{T=0}^{K-1} \alpha_t^2.
\end{align}

Now, in case of decreasing step sizes $\alpha_t=\alpha_0 / \sqrt{t+1}$, we will use the follwing two bounds:

\begin{align}
\sum_{k=1}^K \frac{1}{\sqrt{k}} \geq \sqrt{K}, \quad \sum_{k=1}^K \frac{1}{k} \leq 2+\log K
\end{align}

in the second and third step respectively:

\begin{align}
\min _{0 \leq t<T} \mathbb{E}\left\|\bg_t\right\|_\rho  &
\leq \sum_{t=0}^{T-1} \frac{\alpha_0}{\sqrt{k+1}} \mathbb{E}\left\|\bg_t\right\|_\rho / \sum_{t=0}^{T-1} \frac{\alpha_0}{\sqrt{t+1}} \\
& \leq \frac{1}{\sqrt{T}}\left[\frac{E_f\left(\x^0\right)-E_f^*}{\alpha_0}+\frac{n \bar{K}}{2} \alpha_0 \sum_{T=0}^{K-1} \frac{1}{k+1}\right] \\
& \leq \frac{1}{\sqrt{K}}\left[\frac{E_f\left(\x^0\right)-E_f^*}{\alpha_0}+\alpha_0 n \bar{K}+\frac{\alpha_0 d \bar{K}}{2} \log K\right] \\
& =\frac{1}{\sqrt{K}}\left[\frac{E_f\left(\x^0\right)-E_f^*}{\alpha_0}+\alpha_0 n \bar{K}\right]+\frac{\alpha_0 n \bar{K}}{2} \frac{\log K}{\sqrt{K}}.
\end{align}

In the case of constant step size $\alpha_t=\alpha$, we can simply reformulate the statement to:

\begin{align}
\frac{1}{K} \sum_{T=0}^{K-1} \mathbb{E}\left\|\bg_t\right\|_\rho
\leq \frac{1}{\alpha K}\left[\left(E_f\left(\x^0\right)-E_f^*\right)+\frac{n \bar{K}}{2} \alpha^2 K\right]
=\frac{E_f\left(\x^0\right)-E_f^*}{\alpha K}+\frac{n \bar{K}}{2} \alpha.
\end{align}

§.§ Proof of the fixed point theorem

Before we prove the fixed point theorem, we present the definition of a fixed point as well as reiterating the theorem.
[Fixed point]
$\bs \in \{ \pm 1 \}^n$  is a fixed point for ,
if $\bs^0 = \bs$ in eq. (7) implies that $\bs^t = \bs$ for all $t = 1, 2, ...$, where $\bs^t = \sign(\x^t)$.
We say that  converges if there exists $t < \infty$ such that $\bs^t$ is a fixed point.

Let $\Z^s$ be the normalised Jacobian of the function $f$ at $\mathbf{s} \in \{\pm 1\}^n$ . Let $\x^t$ be the iterative produced by Eq. 7 from the paper. We further assume, that for $i \in [n]$ $(\nabla_{\x}E_f(\x^t))_i \neq 0$ and  $(\nabla_{\x}E_f(\bs))_i \neq 0$. Under the assumption of a diverging learning rate ($\sum_{i \in \mathbb{N}} \alpha_i = \infty$ and $\alpha_i>0$)
, $\mathbf{s}$ is a fixed point for  if and only if ($\Leftrightarrow$) $-\bPi_{\Z^s} (\nabla_{\mathbf{y}}E_f(\mathbf{s})) = \bs$. This point might not exist, in which case we cannot state the convergence of  for any starting point $\x^0 \in \R^n$.

Our proof follows closely the proof of [6]. We first define the fixed point as a stationary point.
We start with the direction $\Rightarrow$. Hence we assume that $\mathbf{s} \in \{ \pm 1\}^n$ is a fixed point.
We start by restating the update rule for point $\x^T$ starting from a fixed point $\bs$:

\begin{align}
\x^T = \x^0 - \sum_{t = 0}^{T - 1} \alpha_i \bPi_{\bZ^t}(\nabla_{\y}E_f(\mathbf{s})). %
\end{align}

Hence $\x^T = \x^0$ holds by definition for all $T\in\mathbb{N}_+$. As the fixed point is the same in the space of reals or binary variables, we are allowed to take the sign of both sides:

\begin{align}
\bs:=\sign(\x^T) = \sign\left(\x^0 - \bPi_{\bZ^s}(\nabla_{\y}E_f(\bs)) \sum_{i = 0}^{T-1} \alpha_i\right).
\end{align}

Taking the limit $T \rightarrow \infty$ and applying the assumption that $\sum_t \alpha_t=\infty$ yields

\begin{align}
s_i=\lim _{T \rightarrow \infty} \operatorname{sign}\left(\x^0-\bPi_{\bZ^s}(\nabla_{\y} E_f(\bs)) \sum_{t=0}^T \alpha_t\right)_i \nonumber
=\sign(-\bPi_{\bZ^s}(\nabla_{\y} E_f(\bs)))_i) =-\bPi_{\bZ^s}(\nabla_{\y} E_f(\bs)))_i.\nonumber
\end{align}

concluding that $\bs=-\bPi_{\bZ^s}(\nabla_{\y} E_f(\bs)))$.

Next, we prove the direction " $\Leftarrow$ ". If $\bs$ obeys that $-\bPi(\nabla_{\y} E_f(\bs))=\bs$ for all $i \in[n]$, then if we take any $\x^0$ such that $\operatorname{sign}\left(\x^0\right)=\bs, \x^t$ will move in a straight line towards the direction of $-\bPi_{\z^s}(\nabla_{\y} E_f(\bs))$. In other words, $\operatorname{sign}\left(\x^t\right)=\operatorname{sign}\left(\x^0\right)=\bs$ for all $t=0,1,2, \ldots$. Therefore, by definition, $\bs$ is a fixed point.

§.§ QUBO term of BM

The binary projection $\bPi_{\bZ^t}(\bv)$ in <ref> admits the following Ising Model or quadratic unconstrained binary optimisation (QUBO) form:

\begin{align}\label{eq:BP_as_qubo}
\bPi_{\bZ^t} (\bv) &= \argmin_{\bg \in \{-1,1 \}^{n}} \bg^\top\sum_{i=1}^m\bQ_i\bg + \bs^\top\bg
\end{align}

where $\bs=-2\sum_{i=1}^m v_i\z_i^t$, $\bQ_i={\z_i^t}{\z_i^t}^\top$ and $\bQ=\sum_{i=1}^m\bQ_i$.

\begin{align}\label{eq:BP}
\bPi_{\bU} (\bv)&= \argmin_{\bg \in \{-1,1 \}^{n}} \sum_{i=1}^m||v_i - \bg^\top \bu_i||_2^2 \\
&= \argmin_{\bg \in \{-1,1 \}^{n}} \sum_{i=1}^m (v_i - {\bu_i}^\top \bg)^\top (v_i - {\bu_i}^\top \bg) \\
&= \argmin_{\bg \in \{-1,1 \}^{n}} \sum_{i=1}^m \bg^\top{\bu_i}{\bu_i}^\top\bg - 2\sum_{i=1}^m v_i {\bu_i}^\top\bg + \sum_{i=1}^m v_i^2 \\
&= \argmin_{\bg \in \{-1,1 \}^{n}} \bg^\top\sum_{i=1}^m\bQ_i\bg + \bs^\top\bg \label{eq:qubo_equation_supp}
\end{align}

where $\bs=-2\sum_{i=1}^m v_i\bu_i$, $\bQ_i={\bu_i}{\bu_i}^\top$ and $\bQ=\sum_{i=1}^m\bQ_i$.

Runtime analysis on quantum and central processing units In general the experiments are relatively easy on the hardware and are run within minutes on a single CPU.
Hence we will only focus on the QPU time, as well as on the graph benchmarks.
For all five Rosenbrock runs we use a total of one minute QPU time.
The small experiments on the Adult dataset as well as the Karateclub network are solved within an hour.
Experiments on graph convolutional neural networks for the bigger datasets Cora and Pubmed run on the CPU for six hours.
In those experiments, we did not need to use any GPU time during training.
The experiments on D-Wave for MLPs are running for 6 hours and one run consumes about a minute of QPU time.
This is explainable with the minor embeddings.
Calculating those minor embeddings is a non-critical combinatorial optimisation problem and, at the moment, is known to be a bottleneck for quantum annealers. We are hopeful that those can be accelerated and make all QA algorithms faster, in the future.

To generate all the experiments in our paper we spent an hour of QPU time on a D-Wave quantum annealer as well as 1000h of CPU time.

§ IMPLEMENTATION DETAILS

Our code, which we will release, is implemented in Pytorch [52].
To calculate all gradients in all the experiments, we use the torch autograd library.
In the graph experiments with Cora and Pubmed we use the Adam optimiser [38] with the pytorch implementation.
For BinaryConnect [21] we use our own implementation. However, contrary to the BinnaryConnect [21] paper we are using our definition of the sign function and not theirs, which projects $\sign(0) = 1$ or does some random assigning.
For ProxQuant we also use our own implementation, following the paper.
The code will be released upon acceptance.

§.§ MNIST Feature Extraction

In the paper, we already described how we extract the features of the MNIST images.
In this section, we give a more detailed overview over this procedure as well as a qualitative example in <ref>.
To this end, we follow a three-step procedure to generate the 16-bit features: (i) we use the key points generated from the method described in [48]. They are marked in <ref> as blue dots in the middle and right image. Afterwards in step (ii) we sample 16 lines, which all run through the centroid of the image. <ref> in the right image contains four out of the 16 lines drawn in green.
For the last and (iii) step, we count the number of feature points above and below the line. If there are more key-points above line $i$ than below, then we assign the feature $i$ to $1$; otherwise to $-1$.

Label generation on MNIST [41] (example).
Left: Initial image of the digit “2”;  Middle: Extracted keypoint features with the superpixel approach [48]; Right: sampling four lines to determine binary features.

§.§ Definitions

Hard tanh
Hard tanh is defined as follows:

\begin{align}\label{eq:hardtanh}
\mathsf{hardTanh}(x) =
\begin{cases}
-1 & \text{ if } x \leq -1, \\
x & \text{ if }  -1 < x < 1, \\
1 & \text{ if } x \geq 1, \\
\end{cases}.
\end{align}

Rosenbrock Function
The Rosenbrock function is defined as follows:

f()=∑_i=1^d-1 f_i(),    where    f_i()=100(x_i+1-x_i^2)^2+(1-x_i)^2

We define the unbiased stochastic gradient of this function as: $\tilde{\bg} = \nabla_{\x} f_i (\x)$, where $i\in \{1, \dots d-1\}$ is uniformly randomly chosen.
To fit this function into our framework, we also define the loss function for this problem.
We use the MSE loss, with a target of 0.
Hence, the minimisation is defined as:

\begin{align}
E(f(\x), 0) = (f(\x) - 0)^2 = f(\x)^2
\end{align}

As the Rosenbrock function is positive the minimum remains the same.

We use the following definition of the sign function:

\begin{align}\label{eq:sign}
\sign(x) =
\begin{cases}
-1 & \text{ if } x < 0, \\
1 & \text{ if } x \geq 0. \\
\end{cases}
\end{align}

§ CDP HYPOTHESIS TEST FOR NEURAL NETWORKS

We will use the $H_0$ hypothesis Z test to determine if our algorithm is more likely to point in the same direction as the gradient.
As a comparison, we will conduct the same test on the established binary baseline signSGD.
We choose as our hypothesis:

$\mathbf{H_0}$ that the projected gradient points in a random direction, that means:

\begin{align}
\begin{split}
\rho_i(\x^t) = \operatorname{Prob}\left( \bPi_{\Z^t}(\tilde{\nabla}_\y E_f(\x^t))_i=\operatorname{sign} (\nabla_\x E_f(\x^t))_i\right) = 0.5
\end{split}
\end{align}

$\mathbf{H_1}$ that the projected gradient points in the same direction

\begin{align}
\begin{split}
\rho_i(\x^t)= \operatorname{Prob}\left( \bPi_{\Z^t}(\tilde{\nabla}_\y E_f(\x^t))_i=\operatorname{sign} (\nabla_\x E_f(\x^t))_i\right) > 0.5
\label{eq:h1hypothesis}
\end{split}
\end{align}

\begin{align}
\begin{split}
\rho_i(\x^t)
= \operatorname{Prob}\left( \bPi_{\Z^t}(\tilde{\nabla}_\y E_f(\x^t))_i=-\operatorname{sign} (\nabla_\x E_f(\x^t))_i\right) > 0.5
\end{split}
\end{align}

We can clearly see that either $\mathbf{H_0}$ or $\mathbf{H_1}$ is true.

The data collection is done by measuring the gradients during training on the MNIST experiment in the paper, where we use the gradient on the full training set as the real non-stochastic gradient.

Let $k$ be the number of occurrences where the gradient points in the same direction as our projected gradient, $p = 0.5$ as our $H_0$ probability and $n$ as the number of compared gradients.
The Z test for binary variables is formulated as follows:

\begin{align}
Z=\frac{k-n p}{\sqrt{n p(1-p)}}.
\label{eq:zvalue}
\end{align}

A Z value in the interval of $[-1.96, 1.96]$ provides a $95\%$ certainty that $\mathbf{H_0}$ is fulfilled.

The following two tables <ref> and <ref> show a sample collection from 5 runs. In the first column, we describe the number of samples, where our projected gradient and the real gradient point in the same direction, in the second column we have all samples and in the last column we note the Z value for the experimental run, which was calculated with <ref>.

As we can see in <ref> and <ref>, when calculating the gradient with  the Z-value is always bigger than $1.96$. This shows that we reject the $\mathbf{H_0}$ hypothesis and accept the $\mathbf{H_1}$ hypothesis.
Hence we can be confident that the CDP holds in our case.
As the tables further present, the projected gradients point in the same direction as the real gradient, then away. Hence we can narrow down the $\mathbf{H_1}$ hypothesis to the CDP assumption.

Further to verify our experiments we also tested signSGD, where the CDP assumption is indicated by a few lemmas. Here we observe that $\mathbf{H_0}$ can be rejected, but on a weaker basis, as the Z values are smaller. This further indicates that our algorithm fulfils the CDP assumption.

$\#$ of same direction     total $\#$ of samples     Z-value

2285     3915     10.4

2199     3647     12.4

2311     3886     11.8

2334     3984     10.8

2356     3990     11.4

Samples from  directions

$\#$ of same direction     total $\#$ of samples     Z-value

2334     4626     0.61

732     1176     8.39

339     630     1.91

2559     4702     6.06

2103     4683     -6.9

Samples from signSGD

§ ADDITIONAL EXPERIMENTS

In this section, we present additional experiments on the Rosenenbrock curve, with continuous weighted MLPs and some D-Wave analysis.

§.§ 2D visualisation of Rosenbrock optimisation

signSGD [58]

Our Q-SPGD (Gurobi)

Our Q-SPGD (D-Wave)

Plotting the path of the optimised inputs of the 2D Rosenbrock function with different optimisers

In this section, we present additionally a qualitative example of the minimisation of the Rosenbrock function.
We compare signSGD with Q-SPGD (Gurobi) and Q-SPGD (D-Wave).
We see in the <ref>, <ref>, <ref>that both optimisers can find the valley of the function, and also follow this valley in the right direction.
However, after 1800 steps signSGD is closer to the global optimum than Q-SPGD (Gurobi).
But we also notice, that  (D-Wave) converges to a similar point like signSGD.
Those results are coherent with the plots in the paper.

§.§ Training of real-valued MLPs

In Def. 2 () we introduced an update rule for continous neural networks.
We will demonstrate now the ability of Q- to train such a continuous MLP.
We employ a simple 2-layered MLP, with an input of 13, hidden and a hidden and output dimension of 3.
As an activation function for the hidden layer we use ReLU and output activation is log softmax.
We train on the whole wine dataset [42] with a rather small batch size of 4.

SGD     signSGD      (Gurobi)

ACC     $0.63 \pm 0.05$     $0.64 \pm 0.05$     $0.67\pm 0.01$

Training accuracy on the Wine dataset over 5 runs. Here mean and stdev

In <ref> we present our mean training accuracy over five runs.
We see that with this setup  fits the neural network better than the strong baselines of SGD and signSGD.
This suggests, that even in the continuous case, our optimiser calculates strong weight updates, which are better suited than the ones calculated by signSGD.
We further point out, that the network on which we trained is relatively small and hence already heavily regularised [65].
We see that this regularisation technique lowers the performance of SGD, but not of .
Therefore we can conclude that our weight updates are also beneficial for optimising continuous neural networks, where only a small neural network can be deployed.

§.§ Additional analysis of the QUBO problems using D-Wave annealer

Left: QUBO problem graph to update the weights of
the second layer in the first training iteration Right: Embedded
problem graph on the Pegasus [12] D-Wave architecture.

To solve QUBO problems on D-Wave, we minor-embed them onto the Pegasus architecture. As shown in <ref>, this operation interprets the logical graph (left) and implements each logical bit with multiple, redundant qubits (right).
The graph in the figure is an exemplary QUBO corresponding to training the third layer in our MNIST experiment.

Evolution of the QUBO
Next, we investigate the development of the QUBO during training. In <ref> we observe how the QUBO changes over time/iterations.
The QUBO matrices shown in the figure are examples from the MNIST experiment (Sec. 3.3.2) taken from the first layer.
We display one QUBO after every $10^{\mathrm{th}}$ update.
We observe here that all QUBOs have a rather densely connected support of the otherwise sparse graph.
This overall sparsity of the graph enables us to embed those with minor embeddings and deploy them to the D-Wave QPU.

Sampling solutions.
After embedding each QUBO problem in the D-Wave Advantage architecture, we query the quantum annealer multiple times and retrieve the solution corresponding to the lowest energy. However, it is still interesting to see the solutions' distribution. To this end, we analyze QUBOs, which calculate the weight updates of the MNIST experiment in the third layer, using D-Wave Explorer and measure a series of QUBOS, which are 10 updates apart. <ref> plots the different energies for the sampled solutions.
The histogram in the upper left corner corresponds to the QUBO update problem in the first iteration and the histogram in the lower right corner corresponds to the QUBO problem after 50 updates (e.g., the $51^{\text{st}}$ update).
In total, we collected 1000 samples from D-Wave and compared the energy.
The resulting distribution can be described as a tail-heavy Poisson distribution, which is shifted towards negative numbers.
This type of resulting distribution is typical for problems, which are difficult to compute.

Similarity In <ref> we compare the similarities of the sampled D-Wave solutions.
Here the 10 scores with the lowest energy are compared with the Jaccard Similarity.
We can see that all QUBO solutions are fairly dissimilar to each other.
As seen in the picture only the first and fourth / fifth solutions have a relatively high similarity score of around $0.4$.
The indications of that are two-fold: (i) it demonstrates the importance of finding the optimal solution, and (ii) even though the D-Wave annealer cannot act like a true Bayesian posterior sampler, the samples from its (modified) posterior are diverse and span the solution space. This also confirms the findings of [11] in which these samples are used as alternative solutions.

Evolution of the QCBO formulation to calculate the weight updates for the first layer with the  algorithm.

Histogram of the different sampling energies of the weight update QUBO of the third layer. In the first column, we show the histogram of the zeroth, tenth, and twentieth updates, and in the second column the histogram of the thirtieth, fortieth, and fiftieth updates from top to bottom
[0.63]Jaccard similarity of the  10 samples with the lowest energy
D-Wave statistics of the energy and similarities of 1000 anneals.

§.§ Ablation studies

Learning rate (lr) comparison.

Batch size comparison.

Impact of learning rate (LR) and batch sizes on the test accuracy. We plot the mean accuracy over six runs.

We now test the effect of different hyperparameters of our algorithm, particularly the batch size and the learning rate.

The effect of batch size
Batching is an important aspect of modern optimisers as it (i) aids theoretical guarantees regarding stochastic convergence, and (ii) becomes essential when training on large datasets. Hence, we now assess the behaviour of our optimiser against different selections of batches. <ref> shows the development of the test accuracy over time as a mean over six runs on the Cora dataset [46]. It is observed that our hybrid optimiser is robust to changes in the batch size, reported as percentages of training data.

Effect of learning rate
Similarly, we plot the mean test accuracy over time for different learning rates in <ref>, which seems to have a bigger impact than the batch size. With a learning rate of either 0.0001 or 0.1, we achieve the best results. However, with learning rates in between, we do not converge to the optimal accuracy within the training duration.

§ EXTENDED RELATED WORK

Binary neural networks (BNNs)
Compared to full-precision neural networks, BNNs [33] consume less memory and provide faster inference at the cost of certain information loss. Much effort has been devoted to mitigating performance degradation due to binarisation, as detailed in recent surveys [55, 75]. Here, we summarise important developments based on the recent surveys. We refer the reader to said surveys for a more complete overview.

Based on their binarisation strategies, BNN models can be categorised into naive or optimisation-based ones. BinaryConnect [21] is an early naive BNN which binarises the weights by a sign function. In this way, most real-valued arithmetic operations are replaced by lightweight bitwise XNOR operations. BinaryNetg [33] is an extension that binarises both weights and activations during inference and training. Dorefa-Net [81] defines new quantisation functions for weights, activations and gradients to accelerate network training.

Optimisation-based BNNs attempt to more directly address the accuracy drop resulting from weight and activation binarisation in the naive models while preserving the compact nature. Following Qin et al. [55], we further categorise such works based on their optimisation goals:

* Minimising the quantisation error. XNOR-Net introduces a scaling factor for the binary weights and activations to generate a better approximation to the corresponding real-valued parameters [56]. XNOR-Net++ builds on this work and fuses the scaling factors of the weights and activations into one. The fused scaling factor is gained from discriminative training via backpropagation [15]. Dorefa-Net defines new quantisation functions for weights, activations and gradients instead of a pre-defined uniform quantiser to accelerate network training [81]. LQ-Nets trains the quantisers for activations in each layer and weights in each channel jointly with the neural network [79]. Parameterised Clipping Activation (PACT) shifts the activation distribution via learning a quantisation scaling factor during network training to achieve low bit-width precision [18]. LAB2 directly minimises binarisation loss for weights using a proximal Newton algorithm [31].

* Improving the loss function. How-to-Train adds a regularisation term to encourage weight bipolarisation in addition to the task-specific loss term [63]. LAB directly develops a loss function during layer-wise weight binarisation [31]. ReActNet introduces a spatial distributional loss to enforce the BNN output distribution to approximate that of the full-precision model, on top of activation distribution shift and reshaping using new quantisation functions [43].

* Reducing the gradient approximation error. Another stream of work focuses on the nature of quantisation and activation functions to enable backpropagation and enhance forward propagation. Bi-Real Net approximates the derivative of the sign function for activations using a piecewise linear function to overcome the intrinsic discontinuity and proposed a magnitude-aware gradient for weights [44]. HWGQ-Net designs a forward quantised approximation and a backward continuous approximation of the ReLU activation function to address the learning inefficiency of BNNs [16].

Various training techniques including asymptotic optimisation, quantisation, gradient approximation and network structure transformation, are adopted in some of the learning-based models above. Standard techniques are used in optimisation, involving Adam, AdaMax, SGD and RMSprop. For more details please refer to Qin  [55].

Quantum annealing
Recent progress in quantum computer vision [11] has shown that various computer vision algorithms could benefit from QUBO-formulations and quantum annealing [11, 17, 62, 78, 2, 47, 76, 28, 9]. Arthur  [3] has similarly shown the advantage for several machine learning models. A new paradigm was developed recently to harness the computational power of quantum annealers for the training of BNNs [59] as well as for optimising in binary variables [76]. As an early work in this direction, a simple version of BNN was examined. As alluded to above, our work improves the scalability of their approach by adopting a hybrid and incremental training scheme.

Quantum neural networks
As an emerging field, quantum machine learning [10, 61, 14, 71] has shown that the training of linear regression, support vector machines and k-means clustering admit QUBO-like formulations [3], whereas the first neural network variants trained on quantum hardware were Boltzmann Machines. In parallel, quantum deep learning [70, 29, 36], the problem of creating quantum circuits that enhance the operations of neural networks by physically realising them, has emerged to alleviate some of the computational limitations of classical deep learning, thanks to the efficient training algorithms [37, 7]. However, a severe challenge remains to be implementing non-linearities and other non-unitary operations with quantum unitaries [60, 36]. Similarly, quantum convolutional neural networks (Q-CNNs) [20] have been developed with the mindset of implementing the analogue of classical CNN operations via quantum gates. These works are fundamentally different to ours, as they try to realise quantum implementations of machine learning algorithms, whereas we leverage quantum computation to solve the classical problem of neural network training.

Network quantisation
Instead of training BNNs from data, network quantisation techniques aim to reduce the precision (including to binary) of a trained full-precision neural network in a manner that preserves accuracy; see Gholami et al. [30] for a recent survey. We stress that quantisation fundamentally differs in formulating a BNN for execution or training on quantum computers.
In the current work, we are focusing on the most extreme kind of network quantisation: binary neural networks. More general types of quantisation can be formulated as a QUBO in principle [49].

Binary graph neural networks
As opposed to the plethora of methods dealing with the binarisation and quantisation of ordinary neural networks, binary graph neural networks are under-explored. Only a few works exist [68, 5, 67] and most of them operate on the principles of XNOR-Net and its variants, aiming to binarise the multiplication. Still, re-quantisation remains to be a necessary operation to avoid full-precision computations [73]. On a parallel track, developing efficient binary graph operators [35] or quantising graph neural networks [82, 35, 32] as well as GNN-compression are somewhat studied. We refer the reader to a recent survey [69] and find it worth mentioning that none of these procedures are friendly for implementation on quantum hardware.