# Using Findroot with Precompiled Function

I am writing a Monte Carlo algorithm to solve some algebraic equations. To do this I need to evaluate FindRoot many times on the same function. Following from F'x's answer to Expression evaluation inside of FindRoot inside a Compiled Function, I understand that FindRoot will compile the target function each time it is evaluated, and so the fastest way for me to work with FindRoot is to combine the function myself beforehand. This should save FindRoot compiling each time I run it.

I can't however work out the syntax to do this. My toy example code is:

f[x_] := x^x + (4 - x)^(4 - x) - 10
g = Compile[{x}, x^x + (4 - x)^(4 - x) - 10];
FindRoot[f[x], {x, 3.99}, Compiled -> False] // Timing
FindRoot[f[x], {x, 3.99}] // Timing
FindRoot[g[x], {x, 3.99}, Compiled -> False] // Timing
FindRoot[g[x], {x, 3.99}] // Timing


None of the ones with g work; I get CompiledFunction::cfsa, and the g ones are slower so I assume this is inputting the uncompiled function.

I have tried also

FindRoot[g, {x, 3.99}] FindRoot[g, {x, 3.99},Compiled-> False]

Which don't evaluate to anything.

What is the correct syntax for this? Or if I've misunderstood the quoted article, could someone explain what I should be doing?

• Try FindRoot[g, {3.99}]. If a function is supplied, no symbol is needed. If you would like to use Newton's method, you have to supply a (compiled) function Dg also for the Jacobian, e.g. with FindRoot[g, {3.99}, Jacobian -> Dg]. For scalar function one may have to fiddle a bit with braces in Dg such that Dg[x] becomes a 1 x 1 matrix; Mathematica will complain otherwise. – Henrik Schumacher Dec 13 '17 at 17:19
• Why do you need to use FindRoot many times on the same function? Won't it always return the same answer? I assume the function must have some changes upon each invocation of FindRoot. Could you clarify what these changes are? – Carl Woll Dec 13 '17 at 20:57
• somebody should make that an answer. @HenrikSchumacher – george2079 Dec 14 '17 at 12:39
• @xzczd I wondered what Chuck Norris would do if he had to implement FindRoot and just gave it a try. – Henrik Schumacher Jan 8 '18 at 15:31
• I think your solution is absolutely intuitive and makes complete sense now that I see it, but I wasn't going to get there just by thinking about Chuck Norris on my own though – Joe Jan 8 '18 at 16:03

## 2 Answers

Finally, I cast my comment into an answer...

To my own surprise, FindRoot[g, {3.99}] works well. If a function is supplied, no symbol is needed. So the current implementation of FindRoot can decide very early that symbolic computation (and maybe further compilation) is futile and it branches to a completely numerical algorithm.

If you would like to use Newton's method, you have to supply a (compiled) function Dg also for the Jacobian, e.g. with FindRoot[g, {3.99}, Jacobian -> Dg]. For a scalar function one may have to fiddle a bit with braces in Dg such that Dg[x] becomes a 1 x 1 matrix; Mathematica will complain otherwise.

See also this post for details.

The key to the CompiledFunction::cfsa error is the option RuntimeOptions -> {"EvaluateSymbolically" -> False}. Maybe this is what you were expecting:

f[x_] := x^x + (4 - x)^(4 - x) - 10
g = Compile[{x}, x^x + (4 - x)^(4 - x) - 10,
RuntimeOptions -> {"EvaluateSymbolically" -> False}];
FindRoot[f[x], {x, 3.99}, Compiled -> False] // RepeatedTiming
FindRoot[f[x], {x, 3.99}] // RepeatedTiming
FindRoot[g[x], {x, 3.99}, Compiled -> False] // RepeatedTiming
FindRoot[g[x], {x, 3.99}] // RepeatedTiming
(*
{0.00047, {x -> 2.37473}}
{0.00048, {x -> 2.37473}}

{0.00036, {x -> 2.37473}}
{0.00036, {x -> 2.37473}}
*)


Without the option setting, the compiled function evaluates symbolically:

g[x]


CompiledFunction::cfsa: Argument x at position 1 should be a machine-size real number.

(*  -10 + (4 - x)^(4 - x) + x^x  *)


FindRoot evaluates the argument symbolically to be able compute the derivative. And it probably uses the symbolic results and then compiles them (or in the case of g[x], recompiles them). So to get the benefit of precompiling g, you need to do something to prevent FindRoot from evaluating g symbolically. One way is the option RuntimeOptions -> {"EvaluateSymbolically" -> False}; the other is to use _?NumericQ as in h.

• Michael, I think a timing comparison between memoized and unmemoized secant method with RepeatedTiming would be only fair if you erase the memoized values of h in between. – Henrik Schumacher Jan 8 '18 at 17:07
• FindRoot[h[x], {x, 3.99}] // RepeatedTiming evaluates FindRoot[h[x], {x, 3.99}] several times and since h is memoized, g will be called only in the first of these runs, since h will already know what its values are afterwards. But that's not the typical situation "in free hunting ground". So, the Clear[h] should also be embraced by RepeatedTiming. An alternative could be to use random numbers as initial guess for FindRoot... – Henrik Schumacher Jan 8 '18 at 17:29
• @HenrikSchumacher Yep, even though the derivative takes twice as long to evaluate than the function, the overhead of FindRoot dominates the timingin this example. Indeed, the overhead of memoizing is measurable. – Michael E2 Jan 8 '18 at 18:04
• Michael, just to clarify: I did not doubt that the secant method can perform better than Newton's method. My comment is just about the memoization overhead. – Henrik Schumacher Jan 8 '18 at 18:23
• @HenrikSchumacher My take-away is that in univariate cases, it's usually not worth taking an extra step to get FindRoot to use Newton's method. – Michael E2 Jan 8 '18 at 18:31