# Findroot with a precompiled function with parameters

If I want to precompile a function that I intend to put into FindRoot many times, I could do it like this:

f[x_] := x^x + (4 - x)^(4 - x) - 10
g = Compile[{x}, Evaluate[f[x]]];
FindRoot[g, {3.99}]

Then if I want to run FindRoot on this function multiple times it will be faster. Generally FindRoot will compile the input function as the function will be called many times, and now if I want to run FindRoot itself many times (for example with different initial conditions to find different roots of the function) the function is compiled only once.

Now suppose that I have a similar function which depends also on a parameter. I'd like to be able to compile the function only once (and not for every value of a), and use FindRoot for different values of the parameter like so:

f[x_, a_] := x^x + (a - x)^(a - x) - 10
FindRoot[f[x, 4], {x, 3.99}]
FindRoot[f[x, 5], {x, 4.99}]

I can't work out I would do this, can anybody help?

In order to profit from compilation, you have to ensure that f gets really compiled. With can help here. In order to allow for Newton's method, we compile the expression for f and its derivative with respect its first argument into CompiledFunctions cf and cDf that expect vector-valued input.

f[x_, a_] := x^x + (a - x)^(a - x) - 10;
Quiet[Block[{x, a},
cf = With[{code = N[{f[x[[1]], a[[1]]]}]},
Compile[{{x, _Real, 1}, {a, _Real, 1}}, code]
];
cDf = With[{code = N[D[{f[x[[1]], a[[1]]]}, {{x[[1]]}, 1}]]},
Compile[{{x, _Real, 1}, {a, _Real, 1}}, code]
];
]
];

We use wrappers F and DF that only evaluate if vector-valued arguments are supplied in order to hand them over to FindRoot. Note the many braces for the initial value.

F[a_] := With[{ccf = cf}, x \[Function] ccf[x, a]];
DF[a_] := With[{ccDf = cDf}, x \[Function] (ccDf[x, a])];
FindRoot[F[{4.}], {{{3.99}}}, Jacobian -> DF[{4.}]]

{{2.37473}}

• OK great thanks I think this is doing what I want. I'd like to be able to package up the F to depend on a so it works something like FindRoot[F[a],{{{3.99}}}], but I can't work that out? Also I'm not quite sure why you need N or to have lots of braces for the initial value, could you explain that? – Joe Jan 8 '18 at 15:01
• N ensures that all integers that pop up in the expression code (e.g. the 10) are cast to doubles at compile time; otherwise there will be type casts from integers to doubles in the compiled code which slows down things unnecessarily. – Henrik Schumacher Jan 8 '18 at 15:13
• The braces are for convincing Mathematica that it should eat a vector (there is one extra pair of braces though that I have to use as a hack which I cannot explain to myself). – Henrik Schumacher Jan 8 '18 at 15:23
• In the meantime, I worked out how to encapsulate it more nicely. Have a look. – Henrik Schumacher Jan 8 '18 at 15:29
• That looks great thanks very much. I think I will need some more time to digest how this works properly so I might come back with a further question later. – Joe Jan 8 '18 at 15:31

I know the specific question is about using a compiled function in FindRoot, but I would like to point out that sometimes compiled functions are slower than the vectorized/auto-parallelized functions in the MKL (on Intel machines) that underlie basic mathematical functions, such as those found in polynomials, in which the OP expressed interest in a comment.

Below is a vectorized FindRoot for solving 100,000 problems of the type in the question. It takes a little over a half second on average. On my machine it uses all cores fully while computing (I have 4, 8 virtual, and the CPU runs at 400%).

ClearAll[f];
f[x_?VectorQ, a_?VectorQ] := x^x + (a - x)^(a - x) - 20;
df0[x_?VectorQ, a_?VectorQ] = D[x^x + (a - x)^(a - x) - 20, x];
df[x_?VectorQ, a_?VectorQ] := DiagonalMatrix@SparseArray@df0[x, a];
params = RandomReal[{3, 5}, 100000];
FindRoot[f[x, params], {x, 2 params/3}, Jacobian :> df[x, params]]; // RepeatedTiming

(*  {0.539, Null}  *)

By comparison, running FindRoot on the OP's nonparametric example 100,000 times takes over 35 seconds. (Even if I use ParallelDo, it takes 8.5 sec., more than 10 times slower.)

ClearAll[f];
f[x_] := x^x + (4 - x)^(4 - x) - 10
g = Compile[{x}, Evaluate[f[x]]];
Do[FindRoot[g, {3.99}], {100000}]; // AbsoluteTiming

(*  {35.7169, Null}  *)

It is hard to be sure from the comment whether such an approach could be adapted to the OP's actual use-case, but it works quite well on the problem described in this question.

• Wow thanks, this is very impressive! So let me see if I can explain this speed increase. 1) FindRoot now only compiles once for each vector of inputs. 2) the calculations are sent to some special part of the processor which is optimised for numerical mathematical operations via this thing called MKL 3) the calculation is parallelised, which means different inputs from the vector are computed on different cores of the processor. Is that all correct? And did I miss anything? – Joe Jan 9 '18 at 11:12
• I need to be able to input initial conditions which aren't so obviously correlated to the choice of parameters. So from your code of FindRoot[f[x, #], {x, 2/3#}, Jacobian -> df[x, #]] &[params] I try extending this in the most obvious way I can think of to FindRoot[f[x, #1], {x, #2}, Jacobian -> df[x, #1]] &[params, inits], but this didn't work. Could you explain how to do this? – Joe Jan 9 '18 at 11:15
• @Joe I'm pretty sure FindRoot does not compile the function at all. Modern CPUs are optimized for arithmetic on arrays. Access in Mathematica is through what are called packed arrays. A number of libraries support this under the hood (MKL, BLAS, LAPACK), which you can google if you're interested. -- Your code FindRoot[f[x, #1], {x, #2}, Jacobian -> df[x, #1]] &[params, inits] worked for me, but I don't know what you used for inits. – Michael E2 Jan 9 '18 at 14:16

A simple idiom to do this is the following:

Module[{fc},
fc[a_] := fc[a] = (Print["Compile called"]; Compile[{x}, x^x + (a - x)^(a - x) - 10]);
f[x_?NumericQ, a_?NumericQ] := fc[a][x];
]

The trick is that when f is called, it calls fc[a] which compiles the function for a specific a and stores it through Memoization. After this first call, it will use the stored compiled function.

Therefore, compile is called only once for each different a.

## Edit

Obviously, I misunderstood the question.

I want to compile only once for all values of a. This was what I meant in my wording of the original question, but I have edited to be more clear.

However, I don't see why you don't simply use

fc = Compile[{x, a}, x^x + (a - x)^(a - x) - 10];
f[x_?NumericQ, a_?NumericQ] := fc[x, a];

It should be noted that your original code

f[x_] := x^x + (4 - x)^(4 - x) - 10
g = Compile[{x}, f[x]];

doesn't do much good. f is not compiled. Instead, you are making the situation worse by creating a compiled function that calls back to the main kernel. Please use <<CompiledFunctionTools and CompilePrint[g] to see your compiled code. Your expression is not included.

Another thing is that the more general expression with a will become complex for certain points and Mathematica will warn you about this.

• I want to compile only once for all values of a. This was what I meant in my wording of the original question, but I have edited to be more clear. – Joe Jan 8 '18 at 13:27
• See my edit then. – halirutan Jan 8 '18 at 13:37
• Huh? Who downvoted here? People go crazy... – Henrik Schumacher Jan 8 '18 at 14:20
• Oh yes right you are my code was doing nothing. I missed the Evaluate in my question I've edited that in now. I'm not really sure how your f[x_?NumericQ, a_?NumericQ] := fc[x, a] answers my problem, I don't see how to put this into FindRoot. Yes thanks for pointing out about the range of the a parameter, this is something that I'm taking care of. – Joe Jan 8 '18 at 15:07
• @Joe Are you referring to that 50µs difference that is spent on evaluating the g[x] (i.stack.imgur.com/Js0pm.png)? I was under the impression your example is an example. Have you tried what happens if you indeed have a complex expression with thousands of terms? Additionally, I wasn't aware of your other post that clears up the background of your question very much. See for instance this example i.stack.imgur.com/7mK9y.png You should consider visiting the Mathematica Chat every once in a while where things can be discussed better. – halirutan Jan 8 '18 at 16:08

Instead of using FindRoot with a parameter, you could use NDSolve. For your example this would be:

f[x_, a_] := x^x + (a-x)^(a-x) - 10

(* initial value *)
x4 = x /. FindRoot[f[x, 4] == 0, {x, 3.99}]

(* ode *)
eqn = D[f[x[a], a], a] == 0

(* NDSolveValue *)
sol = NDSolveValue[{eqn, x[4] == x4}, x, {a, 2, 10}]

2.37473

(a - x[a])^( a - x[a]) (1 + Log[a - x[a]] (1 - Derivative1[x][a]) - Derivative1[x][a]) + x[a]^x[a] (Derivative1[x][a] + Log[x[a]] Derivative1[x][a]) == 0

NDSolveValue::mxst: Maximum number of 51511 steps reached at the point a == 2.450725468002232.

NDSolveValue::ndsz: At a == 4.258744344333925`, step size is effectively zero; singularity or stiff system suspected.

InterpolatingFunction[{{2.45073, 4.25874}}, <>]

And a plot:

Plot[sol[a], {a, Sequence@@sol["Domain"][[1]]}]

• Hi Carl Woll, thanks for an interesting suggestion. So this is a toy example for my full problem which is a set of n polynomial equations in n variables, with n! solutions. I need efficient streamlined code to deal with the rapid increase in number of solutions, is NDsolve likely to offer any performance gain over FindRoot? My intuition says that adding differential equations into the mix would probably slow it down, but I'm open to trying new suggestions – Joe Jan 8 '18 at 16:00
• @Joe The advantage that NDSolve has over FindRoot is that the a very good approximation to the next root is automatically generated. Using FindRoot requires an initial guess, and it is not so simple to always have a good initial guess. – Carl Woll Jan 8 '18 at 16:30