Machine Learning Golf: Multiplication

Question

I'd like to propose a different kind of golfing challenge to this community:

(Artificial) Neural Networks are very popular machine learning models that can be designed and trained to approximate any given (usually unknown) function. They're often used to solve highly complex problems that we don't know how to solve algorithmically like speech recognition, certain kinds of image classifications, various tasks in autonomous driving systems, ... For a primer on neural networks, consider this excellent Wikipedia article.

As this is the first in what I hope to be a series of machine learning golf challenges, I'd like to keep things as simple as possible:

In the language and framework of your choice, design and train a neural network that, given \$(x_1, x_2)\$ calculates their product \$x_1 \cdot x_2\$ for all integers \$x_1, x_2\$ between (and including) \$-10\$ and \$10\$.

Performance Goal

To qualify, your model may not deviate by more than \$0.5\$ from the correct result on any of those entries.

Rules

Your model

• must be a 'traditional' neural network (a node's value is calculated as a weighted linear combination of some of the nodes in a previous layer followed by an activation function),
• may only use the following standard activation functions:
  1. \$\textrm{linear}(x) = x\$,
  2. \$\textrm{softmax}(\vec{x})_i = \frac{e^{x_i}}{\sum_j e^{x_j}}\$,
  3. \$\textrm{selu}_{\alpha, \beta}(x) = \begin{cases} \beta \cdot x & \text{, if } x > 0 \\ \alpha \cdot \beta (e^x -1 ) & \text{, otherwise} \end{cases}\$,
  4. \$\textrm{softplus}(x) = \ln(e^x+1)\$,
  5. \$\textrm{leaky-relu}_\alpha(x) = \begin{cases} x & \text{, if } x < 0 \\ \alpha \cdot x & \text{, otherwise} \end{cases}\$,
  6. \$\tanh(x)\$,
  7. \$\textrm{sigmoid}(x) = \frac{e^x}{e^x+1}\$,
  8. \$\textrm{hard-sigmoid}(x) = \begin{cases} 0 & \text{, if } x < -2.5 \\ 1 & \text{, if } x > 2.5 \\ 0.2 \cdot x + 0.5 & \text{, otherwise} \end{cases}\$,
  9. \$e^x\$
• must take \$(x_1, x_2)\$ either as a tupel/vector/list/... of integers or floats as its only input,
• return the answer as an integer, float (or a suitable container, e.g. a vector or list, that contains this answer).

Your answer must include (or link to) all code necessary to check your results -- including the trained weights of your model.

Scoring

The neural network with the smallest number of weights (including bias weights) wins.

Enjoy!

Answer 1

21 13 11 9 weights

This is based on the polarization identity of bilinear forms which in the one dimensional real case reduces to the polynomial identity:

$$ x\cdot y = \frac{(x+y)^2 - (x-y)^2}{4}$$

So y1 just computes [x+y, x-y] using a linear transformation, and y3 is just the absolute value of y1 as a preprocessing step for the next one: Then the "hard" part is computing the squares which I will explain below, and after that just computing a difference and scaling which is again a linear operation.

To compute the squares I use an exponential series \$s\$ which should be accurate for all integers \$\{0,1,2,\ldots,20\}\$ within around \$0.5\$. This series is of the form

$$ \text{approx_square}(x) = \sum_{i=0}^2 w_i \exp(0.0001 \cdot i \cdot x)$$

where I just optimized for the weights W2 (\$=(w_i)_i\$). This whole approximation comprises again just two linear transformations with an exponential activation sandwiched in between. This approach results in a maximal deviation of about 0.02.

function p = net(x)
% 9 weights
one = 1; 
mone =-1;
zero = 0;
fourth = 0.25;
W1 = [1e-4, 2e-4];
W2  = [-199400468.100687;99700353.6313757];
b2 = 99700114.4299316;
leaky_relu = @(a,x)max(a*x,x); 


% Linear
y0 = [one, one; one, mone] * x;

% Linear + ReLU
y1 = mone * y0;
y2 = [leaky_relu(zero, y0), leaky_relu(zero, y1)];

% Linear
y3 = y2 * [one; one];

% Linear + exp
y4 = exp(y3 * W1); 

% Linear + Bias
y5 =  y4 * W2 + b2;

% Linear
y6 = [one, mone]*y5;
p = y6 * fourth;

end

Try it online!

Answer 2

7 weights

eps = 1e-6
c = 1 / (2 * eps * eps)

def f(A, B):
	e_s = exp(eps * A + eps * B)  # 2 weights, exp activation
	e_d = exp(eps * A - eps * B)  # 2 weights, exp activation
	return c * e_s + (-c) * e_d + (-1 / eps) * B  # 3 weights, linear activation

Try it online!

Uses the following approximate equality for small \$\epsilon\$ based on the Taylor expansion \$ e^x \approx 1 + x + \frac{x^2}{2}\$:

$$ AB \approx \frac{e^{\epsilon A+\epsilon B} - e^{\epsilon A-\epsilon B}}{2 \epsilon^2} - \frac{B}{\epsilon} $$

Picking \$\epsilon\$ small enough gets us within the required error bounds. Note that eps and c are constant weights in the code.

Answer 3

33 31 weights

# Activation functions
sub hard { $_[0] < -2.5 ? 0 : $_[0] > 2.5 ? 1 : 0.2 * $_[0] + 0.5 }
sub linear { $_[0] }

# Layer 0
sub inputA() { $a }
sub inputB() { $b }

# Layer 1
sub a15() { hard(5*inputA) }

# Layer 2
sub a8()  { hard(-5*inputA + 75*a15 - 37.5) }

# Layer 3
sub aa()  { linear(-5*inputA + 75*a15 - 40*a8) }

# Layer 4
sub a4()  { hard(aa - 17.5) }

# Layer 5
sub a2()  { hard(aa - 20*a4 - 7.5) }

# Layer 6
sub a1()  { linear(0.2*aa - 4*a4 - 2*a2) }

# Layer 7
sub b15() { hard(0.25*inputB - 5*a15) }
sub b8()  { hard(0.25*inputB - 5*a8) }
sub b4()  { hard(0.25*inputB - 5*a4) }
sub b2()  { hard(0.25*inputB - 5*a2) }
sub b1()  { hard(0.25*inputB - 5*a1) }

# Layer 8
sub output() { linear(-300*b15 + 160*b8 + 80*b4 + 40*b2 + 20*b1 - 10*inputA) }

# Test
for $a (-10..10) {
        for $b (-10..10) {
                die if abs($a * $b - output) >= 0.5;
        }
}

print "All OK";

Try it online!

This does long multiplication in (sorta) binary, and thus returns the exact result. It should be possible to take advantage of the 0.5 error window to golf this some more, but I'm not sure how.

Layers 1 to 6 decompose the first input in 5 "bits". For golfing reasons, we don't use actual binary. The most significant "bit" has weight -15 instead of 16, and when the input is 0, all the "bits" are 0.5 (which still works out fine, since it preserves the identity inputA = -15*a15 + 8*a8 + 4*a4 + 2*a2 + 1*a1).

Answer 4

43 weights

The two solutions posted so far have been very clever but their approaches likely won't work for more traditional tasks in machine learning (like OCR). Hence I'd like to submit a 'generic' (no clever tricks) solution to this task that hopefully inspires other people to improve on it and get sucked into the world of machine learning:

My model is a very simple neural network with 2 hidden layers built in TensorFlow 2.0 (but any other framework would work as well):

model = tf.keras.models.Sequential([
tf.keras.layers.Dense(6, activation='tanh', input_shape=(2,)),
tf.keras.layers.Dense(3, activation='tanh'),
tf.keras.layers.Dense(1, activation='linear')
])

As you can see, all layers are dense (which most certainly isn't optimal), the activation function is tanh (which might actually be okay for this task), except for the output layer that, due to the nature of this task, has a linear activation function.

There are 43 weights:

• \$(2+1) \cdot 6 = 18\$ between the input and first hidden layer,
• \$(6+1) \cdot 3 = 21\$ between the hidden layers and
• \$(3+1) \cdot 1 = 4\$ connecting the last hidden and the output layer.

The weights have been trained (with an adam optimizer) by a layered fitting approach: First they've been fitted to minimize the mean squarred error not only on integer multiplication between \$-10\$ and \$10\$ but actually on inputs in a certain neighborhood around these values. This results in much better convergence due to the nature of gradient descent. And it accounted for 400 epochs worth of training on 57,600 training samples each, using a batch size of 32.

Next, I've fine-tuned them -- optimizing for the maximum deviation on any of the integer multiplication tasks. Unfortunately, my notes don't show much fine tuning I ended up doing, but it was very minor. In the neighborhood of 100 epochs on those 441 training samples, with a batch size of 441.

These are the weights I ended up with:

[<tf.Variable 'dense/kernel:0' shape=(2, 6) dtype=float32, numpy=
 array([[ 0.10697944,  0.05394982,  0.05479664, -0.04538541,  0.05369904,
         -0.0728976 ],
        [ 0.10571832,  0.05576797, -0.04670485, -0.04466859, -0.05855528,
         -0.07390639]], dtype=float32)>,
 <tf.Variable 'dense/bias:0' shape=(6,) dtype=float32, numpy=
 array([-3.4242163, -0.8875816, -1.7694025, -1.9409281,  1.7825342,
         1.1364107], dtype=float32)>,
 <tf.Variable 'dense_1/kernel:0' shape=(6, 3) dtype=float32, numpy=
 array([[-3.0665843 ,  0.64912266,  3.7107112 ],
        [ 0.4914808 ,  2.1569328 ,  0.65417236],
        [ 3.461693  ,  1.2072319 , -4.181983  ],
        [-2.8746269 , -4.9959164 ,  4.505049  ],
        [-2.920127  , -0.0665407 ,  4.1409926 ],
        [ 1.3777553 , -3.3750365 , -0.10507642]], dtype=float32)>,
 <tf.Variable 'dense_1/bias:0' shape=(3,) dtype=float32, numpy=array([-1.376577  ,  2.8885336 ,  0.19852689], dtype=float32)>,
 <tf.Variable 'dense_2/kernel:0' shape=(3, 1) dtype=float32, numpy=
 array([[-78.7569  ],
        [-23.602606],
        [ 84.29587 ]], dtype=float32)>,
 <tf.Variable 'dense_2/bias:0' shape=(1,) dtype=float32, numpy=array([8.521169], dtype=float32)>]

which barely met the stated performance goal. The maximal deviation ended up being \$0.44350433\$ as witnessd by \$9 \cdot 10 = 90.443504\$.

My model can be found here and you can also Try it online! in a Google Colab environment.

Answer 5

2 weights

I was inspired by the other answers to approximate the polarization identity in a different way. For every small \$\epsilon>0\$, it holds that

$$ xy \approx \frac{e^{\epsilon x+\epsilon y}+e^{-\epsilon x-\epsilon y}-e^{\epsilon x-\epsilon y}-e^{-\epsilon x+\epsilon y}}{4\epsilon^2}.$$

It suffices to take \$\epsilon=0.01\$ for this challenge.

The obvious neural net implementation of this approximation takes weights in \$\{\pm\epsilon,\pm(4\epsilon^2)^{-1}\}\$. These four weights can be golfed down to three \$\{\pm\epsilon,(4\epsilon^3)^{-1}\}\$ by factoring \$\pm(4\epsilon^2)^{-1}=\pm\epsilon\cdot(4\epsilon^3)^{-1}\$. As I mentioned in a comment above, every neural net with weights in machine precision can be golfed to a (huge!) neural net with only two distinct weights. I applied this procedure to write the following MATLAB code:

function z=approxmultgolfed(x,y)

w1 = 0.1;   % first weight
w2 = -w1;   % second weight

k  = 250000;
v1 = w1*ones(k,1);
v2 = w2*ones(k,1);

L1 = w1*eye(2);
L2 = [ w1 w1; w2 w2; w1 w2; w2 w1 ];
L3 = [ v1 v1 v2 v2 ];
L4 = v1';

z = L4 * L3 * exp( L2 * L1 * [ x; y ] );

All told, this neural net consists of 1,250,010 weights, all of which reside in \$\{\pm0.1\}\$.

How to get away with just 1 weight (!)

It turns out you can simulate any neural net that has weights in \$\{\pm0.1\}\$ with a larger neural net that has only one weight, namely, \$-0.1\$. Indeed, multiplication by \$0.1\$ can be implemented as

$$ 0.1x = w^\top wx, $$

where \$w\$ is the column vector of \$10\$ entries, all equal to \$-0.1\$. For neural nets in which half of the weights are positive, this transformation produces a neural net that is \$10.5\$ times larger.

The obvious generalization of this procedure will transform any neural net with weights in \$\{\pm 10^{-k}\}\$ into a larger neural net with the single weight \$-10^{-k}\$. Combined with the procedure in my comment above, it therefore holds that every neural net with machine-precision weights can be transformed into a single-weight neural net.

(Perhaps we should modify how re-used weights are scored in future neural net golfing challenges.)