Polynomial approximation of max function

Question

Let me just say upfront I'm not a mathematician, I'm rather looking for a practical answer to my question. I was wondering if there is a polynomial approximation for the function

$$\max(0,x)=\left\{\begin{array}{ll} x, x>0\\0, x \leq 0\end{array} \right.$$ I was thinking something like the sigmoid function could be useful and since it can be tweaked to output values in the range $[0,1)$, then the answer to my question would be the product $x * sigmoid(x)$, correct? If so, in order to increase accuracy, I would need a sigmoid function with steeper slope around 0, so something like $\frac{1}{1+e^{-kx}}$ for some $k\geq1$. The sigmoid function has then a well-known Taylor series approximation which I could compute in Mathematica.

Answer 1

Limit[x LogisticSigmoid[k x], k -> Infinity] gives ConditionalExpression[x, x > 0] so give a large k and you'll get something close to your Max / Ramp function. Unfortunately a Taylor Series will have poor convergence.

You could use Chebyshev polynomials like this (where I'm using 10 terms, although you can increase this):

f[x_] := Ramp[x]
basis = Array[ChebyshevT[#, x] &, 10, 0];
expansion = 2/Pi Integrate[f[x] #/Sqrt[1 - x^2], {x, -1, 1}] & /@ basis;
approx = -expansion[[1]]/2 + expansion . basis;
Plot[{f[x], approx}, {x, -1, 1}]

HornerForm[approx] gives a fairly compact expression: $$x \left(x \left(\left(x^2 \left(\frac{448}{45 \pi }-\frac{256 x^2}{63 \pi }\right)-\frac{80}{9 \pi }\right) x^2+\frac{40}{9 \pi }\right)+\frac{1}{2}\right)+\frac{1}{9 \pi }$$

Unfortunately, Chebyshev polynomials are only good approximations on the interval [-1, +1].

Answer 2

In addition to @flint's answer, one can "play dumb" and do something like this, finding the least squares difference fit over a range:

With[{f = Ramp, range = {-5, 5}, n = 10},
  With[{params = Array[c, n]},
   With[{eqn = FromDigits[params, x]},
    eqn /. Last@Minimize[
       Integrate[(f[x] - eqn)^2,
        Prepend[range, x]], params]]]] // Together

$$\frac{-21879 x^8+1261260 x^6-26276250 x^4+303187500 x^2+512000000 x+172265625}{1024000000}$$

Plot[{Ramp[x], %}, {x, -5, 5}]

Answer 3

Here's my way to play dumb: Chebyshev interpolation. It's almost equivalent to @flinty's result, but in a different form. (The two approach each other as the degree increases.)

{a, b} = {-1, 2};  (* Pick interval *)
deg = 10;          (* Pick degree *)
xx =               (* Chebyshev nodes *)
  Rescale[Sin[Pi/2 Range[-1. deg, deg, 2]/deg], {-1, 1}, {a, b}];
yy = Max[0, #] & /@ xx; (* y = f(x) *)

(* Barycentric is particularly appropriate and efficient here *)
approx = Statistics`Library`BarycentricInterpolation[
   xx, yy,
   "Weights" -> (* optional: increases accuracy slightly *)
    ReplacePart[
     Table[(-1)^k, {k, 0, Length@xx - 1}], {1 -> 1/2, -1 -> 1/2}]];

Plot[{Max[0, x], approx[x]}, {x, -1, 2}]

You can evaluate the derivative (first derivative only) by passing a 1 as a second argument:

Plot[approx[x, 1], {x, -1, 2}]