How to speed up optimization?

Question

Finding a quadratic polynomial of the best approximation in $L$-metrics to $\sin(x)$ on $[0,5]$, I execute

f[a_?NumericQ, b_?NumericQ, c_?NumericQ] := NIntegrate[RealAbs[Sin[x] - a*x^2 - b*x - c], {x, 0, 5}]; 
NMinimize[{Evaluate[f[a, b, c]],   RealAbs[a] <= 20 && RealAbs[b] <= 20 && RealAbs[c] <= 20}, {a, b, c},
AccuracyGoal -> 5, PrecisionGoal -> 5] // Timing

{3923.94, {0.945758, {a -> -0.22899, b -> 0.700517, c -> 0.278367}}}

As we see

Plot[{Sin[x], a*x^2 + b*x + c /. {a -> -0.22898986824771952, 
b -> 0.7005170433410821, c -> 0.27836651693875075}}, {x, 0, 5}]

, the result is satisfactory, but it takes a lot of time. How to speed up that execution and similar ones (say, with dozen variables)? I don't find out an answer here.

Answer 1

This is a bit out there, but you can actually do this using the neural network facilities which opens up the possibility to use a GPU:

Remove["Global`*"];
f = Sin[#] &; (* Our target )
g = {1, #, #^2} &; ( various functions to sum with learned weights to approximate f. *)
net = NetGraph[{
    FunctionLayer[g, "Input" -> "Real", "Output" -> Length[g[0]]],
    NetArrayLayer["Output" -> Length[g[0]]],
    DotLayer[],
    FunctionLayer[-f[#] &],
    TotalLayer[],
    FunctionLayer[RealAbs]}, {NetPort["Input1"] -> 1, 1 -> 3, 2 -> 3, 
    NetPort["Input1"] -> 4, 4 -> 5, 3 -> 5, 5 -> 6}];
net = NetInitialize[net];
training = # -> 0 & /@ RandomReal[{0, 5}, 100000];
net = NetTrain[net, training, TargetDevice -> "GPU", BatchSize -> 256];
coeffs = NetExtract[net, 2][];
Plot[{f[x], coeffs . g[x]}, {x, 0, 5}]

The resulting coefficients coeffs are the optimized {c,b,a} in your quadratic, so the polynomial is given by the dot product coeffs . {1, x, x^2} You may wonder where the integration has gone. Well I believe the gradients are accumulated over the whole batch, so the overall change to the parameters is based on evaluation at many random points in the domain, not just one at a time. In this way, the effect of batching is similar to Monte Carlo integrating the loss first, then updating the parameters so they produce a lower loss for the overall batch.

---

You could also consider not doing the integration, but instead sampling at random points. This cuts down the time by a lot, and decreases the accuracy. It is somewhat similar to the neural method above, but uses global optimization instead and only requires a small adjustment to your code:

f[a_?NumericQ, b_?NumericQ, c_?NumericQ] :=
 Sum[RealAbs[Sin[x] - a*x^2 - b*x - c], {x, RandomReal[{0, 5}, 500]}, 
  Method -> "Procedural"]
{err, sol} = 
 Quiet[NMinimize[{Evaluate[f[a, b, c]], 
    RealAbs[a] <= 20 && RealAbs[b] <= 20 && RealAbs[c] <= 20}, {a, b, 
    c}, MaxIterations -> 100, 
   Method -> {"DifferentialEvolution", "PostProcess" -> False}],
  NMinimize::cvmit]

Extending to other norms is easy too, just define f like this, where p is a number, or $\infty$ for the $L_\infty$ max norm:

f[a_?NumericQ, b_?NumericQ, c_?NumericQ] :=
 Norm[Sin[#] - a*#^2 - b*# - c & /@ 
   RandomReal[{0, 5}, 500], p]

Answer 2

Have you considered using Orthogonalize for this? This finds an orthogonal basis for the powers of x over [0,5] and then projects Sin onto it to find the coefficients in that basis. The result looks the same and it's nearly instantaneous for the quadratic case:

innerProduct[f_, g_] := Integrate[f g, {x, 0, 5}]
basis = Orthogonalize[Table[x^k, {k, 0, 2}], innerProduct];
result = 
 Sum[Projection[Sin[x], basis[[i]], innerProduct], {i, 
   Length[basis]}]
Plot[{result, Sin[x]}, {x, 0, 5}]

You get exact coefficients:

CoefficientList[result, x]
(* {2/5 Sin[5/2]^2 + 3/25 (5 + 5 Cos[5] - 2 Sin[5]) + 
  1/25 (13 - 13 Cos[5] + 30 Sin[5]), -(6/
    125) (5 + 5 Cos[5] - 2 Sin[5]) - 
  6/125 (13 - 13 Cos[5] + 30 Sin[5]), 
 6/625 (13 - 13 Cos[5] + 30 Sin[5])} *)
(* {0.365393, 0.533721, -0.186771} *)

---

This is $L_2$ only though, and unfortunately I don't think you can define an inner product with an $L_1$ norm as it would violate the parallelogram law.