# Polynomial approximations of nonlinearities in neural networks

Imagine, that the only operations I have are scalar addition and scalar multiplication and I want to implement different nonlinearities for neural networks with them. The only option I see here is to use some polynomial approximations.

I can use Taylor series approximations, but it gives poor results. For example, for smooth ReLU $\mathrm{ln}(1 + e^x)$ its Taylor approximation with 9 terms greatly diverges outside the range $[-4,4]$ (look at the plots on WolframAlpha). This is unacceptable because activations in neural networks can go far beyond this range. For other nonlinearities (logistic sigmoid and tanh) we have a similar picture. Of course, I can use much more terms, but it will be computationally intractable.

So, I have two questions:

1. Are there any good polynomial approximations of nonlinearities used in neural networks (ReLU, logistic, tanh)?
2. What range should I care about? For example, if some function $f(x)$ is a good approximation of ReLU in the range [-100, 100], will it be enough for me in practice? In what range are activations of nonlinearities in AlexNet, GoogLeNet, VGG, etc?
• It seems that in general the range of activations is entirely data-dependent. A researcher could undertake to constrain the range of activations in some way, but if that's the case, the range of activations is already known to the researcher because s/he devised the range!
– Sycorax
Aug 22, 2017 at 16:34
• @Sycorax: The idea of Batch Normalization, along with proper weight initialization (e.g. Xavier) largely solves this. You're dynamically shifting/scaling each layer's output, against an accruing mean and standard deviation, along with providing an initialization that semi-intelligently anticipates the output of each layer. May 22, 2018 at 20:48
• @Wunsch Punsch: It's not really clear what your motivation here is. Are you trying to train a network from scratch using polynomial activations, or are you trying to approximate an existing (already trained) network with polynomial activations? If it's the former, then proper layer normalization will somewhat help your neurons, although when dealing with cubics and upward, you run into the squashing/blow up issue. So perhaps you could use the theory of gradient clipping and reflection (and rmsprop techniques). May 22, 2018 at 20:50
• @AlexR. I wanted to compute all the network in TT-format. One can compute matrix-matrix multiplications and additions in TT-format. But to compute the activations we have to unpack the whole thing and leave the format, which is expensive. After computing the activations we pack our network again into it. So polynomial activations could help to avoid this. I think, that BN won't help, because our preactivation values will be spread on a large range anyway and so Taylor approximation will be quite inaccurate May 29, 2018 at 17:36
• So to be crystal clear, are you already given a neural network with weights and activations, and you now want to approximate each activation with polynomials, to then compute it in TT format? May 29, 2018 at 18:15

The problem you're having is due to the asymptotic behavior of the remainder between the Taylor approximation of a function, and the function itself.

If $f$ is at least $k$ times differentiable and you approximate it with a $k$'th order polynomial, $$P_k(x) = f(a) + f'(a) (x - a) + \frac{1}{2} f''(a) (x-a)^2 + \dots + \frac{1}{k!} f^{(k)}(a) (x-a)^k,$$

Then the scaling of the error $P_k(x) - f(x)$ goes as $o(|x - a|^k).$ I'm using the little-o notation. This is particularly noticeable for an infinite-differentiable function, such as the smooth-RELU.

The example that's often used to illustrate this concept is $\sin x.$ The function $\sin x$ has infinite stationary points. This is required for the constant oscillation behavior. However, a polynomial of order $k$ can only have, at max, $k-1$ stationary points. So trying to fit a polynomial to a sine function can approximate it decently well around the point you're referencing. In fact, the higher order your polynomial, the wider the domain of validity your approximation will be. But it will always eventually diverge. The example you gave, the smooth RELU, also has its own assortment of problems. You can either simulate it well near the origin, or at some infinite asymptotic point, but not both. To see this, look at its asymptotic behavior as $x \rightarrow \infty:$

$$\begin{split} \lim_{x \rightarrow \infty} \log (1 + e^x) &= \lim_{x \rightarrow \infty} \log \left[\frac{e^{-x} + 1}{e^{-x}}\right] \\ &= \lim_{x \rightarrow \infty} \log \left[ \frac{1}{e^{-x}} \right] \rightarrow x. \end{split}$$ That is, its infinite asymptotic behavior is just that of the function $g(x) = x.$ However, there is an inflection point near the origin, which cannot possibly be represented well by a first order polynomial, so you need a higher order polynomial to capture that. However, if you use a polynomial $P_k(x)$ of order $k>1$ then, $$\lim_{x \rightarrow \infty} \frac{P_k(x)}{x} = \infty.$$ That is, you cannot possibly have consistent behavior as you increase $x.$ The higher order polynomial will always overtake $x$ at some point, and continue to grow more rapidly, whereas $x$ is the asymptotic behavior you're looking for. Furthermore, its asymptotic behavior on the left side is $h(x) = 0,$ which is itself a polynomial. So you have zero'th order behavior on the left, first order behavior on the right, and an inflection point near the middle, which suggests a polynomial of at least order 2. You cannot possible fit all of these behaviors by a single polynomial.

It is a matter of scale. Looking at the sine wave approximation, the polynomial is indeed a solid approximation of the function, and in fact, this is used by most calculators to give you the output of a $\sin$ function. But this only works within a certain range. If you want to approximate activation functions in your Neural Network you have to keep in mind that the input to the activation functions are somewhat uncontrolled, because they're the linear combination of input nodes that you're fitting by looking for the appropriate weights. You have to ensure that your weights are small enough that the domain of validity of your Taylor approximations hold. Regularization, such as $L2$ weight decay, might do the trick, but it's no guarantee. Another way around this is to increase the order of the polynomial, to make the domain of validity larger. But, given the asymptotic arguments I illustrated above, there is no general way to ensure that a polynomial gives you the behavior you expect out of an infinitely-differentiable activation function.

So the issue is that Taylor series are not always the best finite approximation for a function. With the smooth relu example, you can take advantage of a few tricks to make the approximation better. As an example, $\ln(1+e^{x})=\ln(e^x(e^{-x}+1))=x+\ln(1+e^{-x})$. This form is advantageous when $x$ is large, because then $\ln(1+e^{-x})=e^{-x}-e^{-2x}/2+...$ where all the next terms drop off exponentially.

If you're willing to divide by polynomials, then an even better series is:

$\ln(x) = 2\sum_{k=0}^\infty\frac{1}{2k+1}\left(\frac{x-1}{x+1}\right)^{2k+1},$

which usually gives you great approximations with only a few terms, as long as $x$ is near 1 (but say, x=2 isn't terrible either). By using $\ln(xy)=\ln(x)+\ln(y)$ and $\ln(x)=a\ln(x^{1/a})$, you can always reduce $x$ to something near 1 and then apply the above approximation. By using batch normalization, you can guesstimate the variance of your inputs and pick an $a$ that does the trick for most inputs.