Data mining and Machine Learning research has led to a wide variety of training methods and algorithms for di erent types of models. Many of these methods solve or approximate NP-hard optimization problems at their core, using vastly di erent approaches, some algebraic, others heuristic. This paper demonstrates another way of solving these problems by reducing them to quadratic polynomial optimization problems on binary variables. This class of parametric optimization problems is well-researched and powerful, and o ers a unifying framework for many relevant ML problems that can all be tackled with one e cient solver. Because of the friendly domain of binary values, such a solver lends itself particularly well to hardware acceleration, as we further demonstrate in this paper by evaluating our problem reductions using FPGAs.
Hardware acceleration for machine learning usually involves GPU implementations that can do fast linear algebra to enhance the speed of numerical computations. Our approach is di erent to GPU programming in that we make use of a xed class C of parametric optimization problems that we solve directly on e cient specialized hardware. Solving a problem instance in C thus amounts to nding the correct parameters and feeding them to the solver.
The underlying idea of using a non-universal compute-architecture for
machine learning is indeed not new: State-of-the-art quantum annealers rely on the
very same principle of parametric problem classes. There, optimization problems
are encoded as potential energy between qubits; the global minimum of a loss
function can be interpreted as the quantum state of lowest energy [
Nevertheless, the particular class of optimization problems that quantum annealers can solve is well understood which motivates its use for hardware accelerators outside of the quantum world. 2
A pseudo-Boolean function (PBF) is any function f : Bn 7! R that assigns a
real value to a xed-length binary vector. Every PBF on n binary variables can
be uniquely expressed as a polynomial of some degree d n with real-valued
coe cients [
It has been shown that all higher-degree pseudo-Boolean optimization
problems can be reduced to quadratic problems [
If no specialized algorithm is known for a particular hard combinatorial optimization problem, randomized search heuristics, like simulated annealing or evolutionary algorithms (EA), provide a generic way to generate good solutions.
Inspired by biological evolution, EAs employ recombination or mutation on a
set of \parent" solutions to produce a set of \o spring" solutions. A loss function,
also called tness function in the EA-context, is used to select those solutions
which will constitute the next parent generation. This process is repeated until
convergence or a pre-speci ed time-budget is exhausted [
Motivated by the inherently parallel nature of digital circuits, we developed a highly customizable ( + )-EA FPGA hardware architecture implemented in VHDL. Here, customizable implies that di erent types and sizes of FPGA hardware can be used.
Moreover, our hardware synthesizer allows the end-user to customize the maximal problem dimension n, the number of parent solutions , the number of o spring solutions , and the number of bits per coe cient ij . In case of low-budget FPGA, this allows us to either allocate more FPGA resources for parallel computation ( and ) or for the problem size (n and ).
We used this hardware optimizer for evaluating the QUBO and Ising model embeddings we are going to present in the following sections. 4
With an e cient solver for QUBO and the Ising model at hand, solving a speci c optimization problem reduces to nding an e cient embedding into the domain of n-dimensional binary vectors and devising a method for calculating such that the minimizing vector x corresponds to an optimal solution of the original problem. The actual optimization step then amounts to uploading to the hardware solver which approximates a global optimum that yields an approximately optimal solution to the original problem, once decoded back into the original domain. As we will show, it is possible for multiple valid embeddings to exist for the same problem class, but their e ectiveness di ers from instance to instance. 4.1
A prototypical data mining problem is k-means clustering which is already
NPhard for k = 2. To derive the coe cients for an Ising model solving the
2-means clustering problem, we can use the method devised in [
First, we assume that D has zero mean, such that Pi xi = 0. Next, we compute the Gramian matrix G, where every entry corresponds to an inner product of two data points:
Gij = hxi; xj i
As shown in [
For simplicity we will stick to the linear kernel and assume that G is already centered. As our hardware optimizer encodes as signed b-bit integers (with b 32), we have to round the entries of G, which are real-valued. To minimize loss of precision, we scale the parameters to use the full range of b bits before rounding. As the overall objective function is a linear combination of , scaling all coe cients by is the same as multiplying the objective function by , which means that the position of the optimum is una ected. The nal formula for a single coe cient comes out to ij = b Gij + 0:5c with = Notice that there is no practical di erence in precision between integer and xed-point representation, as the latter is nothing more than a scaled integer representation to begin with.
Exemplary optimization results using this method on the UCI datasets Iris and Sonar are shown in Fig. 1.
20 40 60 80 100 20 40 60 80 100
Time (ms) 200 400
600
Time (ms)
800
1000
200
400
600
800
1000
A linear Support Vector Machine (SVM) is a classi er that { in the simplest
case { takes a labeled dataset D = f(xi; yi) j 1 i ng X Y of two classes
(Y = f+1; 1g) and tries to separate them with a hyperplane [
As for real-world data perfect linear separability is unlikely, each data point is assigned a slack variable i such that i > 0 indicates that xi violates the correctness condition. This yields the second objective, which is to minimize the total slack, so that the entire minimization problem comes out to be minimize 12 kwk22 + C X i s.t. 8i: (hw; xii i b) yi 1 i
The optimization is done by adjusting w, b and i, while C is a free parameter controlling the impact of wrong classi cation on the total loss value. In fact, C is nothing more than an inverse regularization factor, as the loss function can be rewritten as Pi i + kwk22, which illustrates that SVM learning is merely a regularized hinge loss minimization. Typically this problem is solved using its Lagrange dual, which is
n maximize X
i i=1 s.t. 8i: 0
n n 1 X X 2
Note that eq. (3) has striking similarity to eq. (1), the QUBO objective function, only a) it is a maximization instead of a minimization problem, b) index j goes from 1 to n instead to i, which makes for a full matrix instead of a lower triangle matrix, and c) i is continuous between 0 and C instead of a boolean. Points a) and b) are easy to deal with by ipping the sign and multiplying all entries except for the main diagonal of the triangle matrix by 2. Point c) requires a radical simpli cation: Instead of a continuous i 2 [0; C] we assume that i is either 0 or C, reducing it to a boolean indicator whether xi is a support vector or not. Thus we can substitute i for C zi with zi 2 B and end up with the following QUBO formulation: n X i=1
n 0i 1 1
Czi + Xi=1 @Xj=1 C2tij zizj + 12 C2tiiziA where tij = yiyj hxi; xj i In terms of eq. (2) we can express ij = ( 1 C2tii 2 C2tij as C if i = j otherwise
Another typical NP-hard ML problem is to determine the most likely con
guration of variables in a Markov Random Field, known as the MAP prediction
problem [
For the rst QUBO embedding of the MAP problem we encode the assignments of all Xv as a concatenation of one-hot vectors (X1 = xi11 ; : : : ; Xm = ximm ) 7! 0 : : : 010 : : : 0 : : : 0 : : : 010 : : : 0; | {z } | {z } jX1j jXmj where m = jV j and xik is the i-th value in Xk. The QUBO problem dimension is equal to the total number of variable states, n = Pv2V jXvj. The weights are encoded into the quadratic coe cients: Let : V X 7! N be a helper function that, given a vertex and variable state, returns the corresponding bit index within the concatenated one-hot vector. Given two indices i = (u; x) and j = (v; y) with u 6= v, the corresponding coe cient is ij = uv=xy; the negative sign turns MAP into a minimization problem. If u = v and i 6= j, the indices belong to the one-hot encoding of the same variable; as a variable can only be in one state at a time, these two bits cannot both be 1, as this would lead to an invalid one-hot encoding. Therefore the corresponding coe cient must be set to an arbitrary positive penalty value P , large enough to ensure that even the worst valid encoding has a lower loss value than P . On the other hand, if i = j then ij = 0, because the variable assignments have no weights on their own.
ij = 80 > <
P > : where i = (u; x); j = (v; y) uv=xy if i = j if i 6= j and u = v otherwise
The dimension can be reduced by removing bits with index i where all ij are either 0 or P for all j i, as those can never be part of an optimal solution.
For a di erent QUBO embedding, we may assign bits zuv=xy to all non-zero weights uv=xy between speci c values x and y of two variables Xu and Xv, indicating that Xu = x and Xv = y (see g. 2). Again, to avoid multiple or con icting variable assignments we introduce penalty weights between pairs of edges that share a variable but disagree in its speci c value (see g. 3).
The result z of the optimization is a vector of indicators for \active" edges indexed by the set of all non-zero weights , from which we can derive an assignment for every variable, provided it has at least one incident edge of non-zero weight:
Xv = [fy~ j u 2 V nfvg; u 2 Xu; y~ 2 Xv; zuv=xy~ = 1g For a valid encoding, the above set is always a singleton, such that its union yields the contained element. 5
We evaluated our embeddings using our optimizer hardware on multiple UCI datasets.
For 2-means clustering we took ve datasets and compared the resulting
k-means loss values of the Ising model embedding to Lloyd's algorithm [
The SVM embedding was tested extensively on the Sonar dataset, where we created ten random splits of 2=3 training data and 1=3 test data. Using our hardware optimizer, we then trained the binary SVM on each test set of each split with C = 10k for k 2 f 3; 2; 1; 0; 1; 2g and calculated the accuracy on the respective test set; see table 2 for the full results. For C = 10 3, our optimizer found the optimum z = 1, which means that every data point was a support vector. For bigger C the number of support vectors decreased, until for C = 102 the optimum was z = 0, which is why 102 is missing from the table, as w and b could not be calculated without at least one support vector. The best accuracy was achieved with C = 1 and came out to be about 75%. On the same splits we trained a conventional SVM (libsvm using the svm method from R's e1071 package3) and found that the test accuracies were in fact very similar, the di erence being less than one percent on average (see table 3), which is a promising result considering the radical simpli cation steps taken for this model.
The two MRF embeddings were tested on the Sonar and Mushroom datasets, using weights from pre-trained integer MRFs with Chow-Liu tree structures. Let j j denote the number of non-zero weights of a trained MRF: While the MRF on Sonar had only few non-zero weights (j j = 102), Mushroom yielded a much more densely connected MRF (j j = 679). Convergence plots for Mushroom can be seen in g. 5. Obviously the state embedding (where each bit encodes one variable state as a one-hot vector) is superior on Mushroom, as it converges much more quickly to a very good optimum, whereas the edge embedding (where each bit encodes one value from ) takes much longer and does not converge to 3 https://www.rdocumentation.org/packages/e1071/versions/1.7-1/topics/svm nearly as good an optimum, even after several minutes. On the contrary, the edge embedding on Sonar is much more e cient than the state embedding (see g. 6). The di erent convergence rates are due to the di erent QUBO dimensions of the embeddings: The state embedding on Mushroom has dimension ns = Pv2V jXvj = 127, which can be further reduced to 107 by removing unnecessary bits as described in section 4.3. The edge embedding on the other hand leads to a dimension ne = j j = 679. As our hardware optimizer needed about O(n2) to perform the optimization, the latter encoding lead to a performance about 40 times slower. The Sonar MRF has ns = 662, which the reduction step improves signi cantly to 117. As j j = 102 (and ne = 102, consequently) for the Sonar MRF, the edge embedding is even better, though.
As a general rule we observe that for relatively dense MRFs, the (reduced) state encoding leads to faster convergence, while for sparse MRFs with j j < Pv2V jXv j (where Xv is the set of variable states that have at least one non-zero associated with them) the edge encoding is preferable. 6
QUBO and Ising models constitute simple but versatile optimization problems that capture the essence of many problems highly relevant to Machine Learning, and in this work we have summarized some interesting embeddings of NP-hard Machine Learning problems into both models. A general hardware-based solver for either of these models is thus a powerful tool for Machine Learning, and the domain of bit vectors makes for an e cient problem representation even when resources are restricted, e.g. in embedded systems. Further we showed that there can exist multiple embeddings of the same problem that perform di erently depending on properties of the problem instance.
It would be interesting to uncover more embeddings of ML problems in the future, both with conventional hardware acceleration like we did, but also with the advent of quantum annealing in mind. κ = 1/n κ = 2/n κ = 3/n κ = 4/n κ = 5/n κ = 6/n κ = 7/n κ = 8/n κ = 9/n κ = 10/n 0
Acknowledgement This research has been funded by the Federal Ministry of Education and Research of Germany as part of the competence center for machine learning ML2R (01jS18038A)