Find the largest root of a polynomial with a neural network
14,674,000,667 5,436,050 5,403,448 10,385 5,994 4,447  3,806 total precision

For a baseline, I investigated the following approach: Select \$M,\delta,\epsilon&gt;0\$ such that if we sample the polynomial \$p(x)=x^3+ax^2+bx+c\$ at 

$$ S:=\{-M,-M+\delta,-M+2\delta,\ldots,M\}, $$

then the largest sample point \$s^\star\in S\$ satisfying \$p(s^\star)&lt;\epsilon\$ necessarily exists and necessarily resides within \$0.1\$ of the largest root of \$p\$. One can show that for our collection of polynomials, one may take \$M=11\$, \$\delta=0.1\$, and \$\epsilon=10^{-4}\$. 

To design a neural net that implements this logic, we start with a layer of neurons that sample the polynomial on \$S\$. For each \$s\in S\$, we take

$$ x_{1,s} = s^2\cdot a + s\cdot b + 1\cdot c + s^3. $$

Next, we identify which of these are less than \$\epsilon=10^{-4}\$. It turns out that for \$s\in S\$, it holds that \$p(s)&lt;10^{-4}\$ only if \$p(s)\leq 0\$. As such, we can use relu activations to exactly identify our samples:

$$ \frac{\mathrm{relu}(10^{-4}-t) - \mathrm{relu}(-t)}{10^{-4}} = \left\{\begin{array}{cl}1&amp;\text{if }t\leq 0\\0&amp;\text{if }t\geq 10^{-4}.\end{array}\right.$$

We implement this with a few layers of neurons:

$$ \begin{aligned}
x_{2,s} &amp;= \mathrm{relu}(-1\cdot x_{1,s}+10^{-4}), \\ 
x_{3,s} &amp; = \mathrm{relu}(-1\cdot x_{1,s}), \\ 
x_{4,s} &amp;= 10^{4}\cdot x_{2,s}-10^{4}\cdot x_{3,s}. \end{aligned}$$

At this point, we have \$x_{4,s}=1\$ when \$p(s)&lt;10^{-4}\$, and otherwise \$x_{4,s}=0\$. Recall that we seek the largest \$s^\star\$ for which \$x_{4,s^\star}=1\$. To this end, we label \$x_{4,M}\$ as \$x_{5,M}\$ (for notational convenience), and for each \$k\geq 1\$, we iteratively define

$$ x_{5,M-k\delta} = 1\cdot x_{4,M-k\delta}+2\cdot x_{5,M-(k-1)\delta} = \sum_{j=0}^k 2^{k-j}x_{4,M-j\delta}. $$

Thanks to this transformation, every \$x_{5,s}\$ is a nonnegative integer, and \$s^\star\$ is the unique \$s\$ for which \$x_{5,s}=1\$. We may now identify \$s^\star\$ by another application of relu activations:

$$ \mathrm{relu}(t-2)-2\cdot\mathrm{relu}(t-1)+t = \left\{\begin{array}{cl}1&amp;\text{if }t=1\\0&amp;\text{if }t\in\mathbf{Z}_{\geq 0}\setminus\{1\}.\end{array}\right.$$

Explicitly, we define neurons by

$$ \begin{aligned}
x_{6,s} &amp;= \mathrm{relu}(1\cdot x_{5,s} - 2), \\
x_{7,s} &amp;= \mathrm{relu}(1\cdot x_{5,s} - 1), \\
x_{8,s} &amp;= 1\cdot x_{6,s} - 2\cdot x_{7,s} + 1\cdot x_{5s}.
\end{aligned} $$

Then \$x_{8,s}=1\$ if \$s=s^\star\$ and otherwise \$x_{8,s}=0\$. We linearly combine these to produce our output node: 

$$ x_9 = \sum_{s\in S} s\cdot x_{8,s} = s^\star. $$

For the score, each layer has neurons with different levels of precision:
(1) \$6+3+1+9=19\$, (2) \$1+4=5\$, (3) \$1\$, (4) \$5+5=10\$, (5) \$1+1=2\$, (6) \$1+1=2\$, (7) \$1+1=2\$, (8) \$1+1+1=3\$, (9) \$3|S|\$. Furthermore, all but two of the layers have \$|S|=221\$ neurons; layer 5 has \$|S|-1\$ neurons and layer 9 has \$1\$ neuron.

Edit: Improvements: (1) We can sample the polynomial much more efficiently using finite differences. (2) We can bypass layers 2 through 4 by instead using a sigmoid activation. (3) The overflow issues in layer 5 can be averted (and subsequent layers can be combined) by more carefully applying relu activations. (4) The final sum is cheaper with summation by parts.

What follows is the MATLAB code. To be clear, prec is a function (found here) that computes the precision of a vector of weights or biases.

function sstar = findsstar2(a,b,c)

relu = @(x) x .* (x&gt;0);

totprec = 0;

% x1 samples the polynomial on -11:0.1:11
x1=[];
for s = -11:0.1:11
    if length(x1) &lt; 5
        w1 = [s^2 s 1];
        b1 = s^3;
        x1(end+1,:) = w1 * [a; b; c] + b1;
        totprec = totprec + prec(w1) + prec(b1);
    else
        w1 = [-1 4 -6 4];
        x1(end+1,:) = w1 * x1(end-3:end,:);
        totprec = totprec + prec(w1);
    end
end

% x4 indicates whether the polynomial is nonpositive
w4 = -6e5;
b4 = 60;
x4=[];
for ii=1:length(x1)
    x4(end+1) = sigmf(w4 * x1(ii) + b4, [1,0]);
    totprec = totprec + prec(w4) + prec(b4);
end

% x6 indicates which entries are less than or equal to sstar
x5 = zeros(size(x1));
x6 = zeros(size(x1));
x5(end) = 0;
x6(end) = 0;
for ii = 1:length(x5)-1
    w5 = [-1 -1];
    b5 = 1;
    x5(end-ii) = relu(w5 * [x4(end-ii); x6(end-ii+1)] + b5);
    totprec = totprec + prec(w5) + prec(b5);
    w6 = -1;
    b6 = 1;
    x6(end-ii) = w6 * x5(end-ii) + b6;
    totprec = totprec + prec(w6) + prec(b6);
end

% a linear combination produces sstar
w7 = 0.1*ones(1,length(x1));
w7(1) = -11;
sstar = w7 * x6;

%disp(totprec) % uncomment to display score

end