Expression evaluation inside of FindRoot inside a Compiled Function

Question

I'm trying to get some performance increase out of my own implicit differential equation solver using Compile[]. The uncompiled function is of the following form:

fun = Module[
  {sol = ConstantArray[0, 10], init = 0.001},
  f[xN_] = xN + xNP1^2 == 4;
  Do[
    sol[[i]] = xNP1 /. FindRoot[f[init], {xNP1, init}];
    init = sol[[i]]
    , {i, 10}
  ];
  sol
]

which works correctly. Of course, for the real function I need many more than 10 iteration in the loop and was hoping to gain some performance increase with Compile[]. Here is the compile code:

cFun = Compile[{},
         Module[
           {sol = ConstantArray[0, 10], init = 0.001},
           f[xN_] = xN + xNP1^2 == 4;
           Do[
             sol[[i]] = xNP1 /. FindRoot[f[init], {xNP1, init}];
             init = sol[[i]]
             , {i, 10}
           ];
           sol
         ]
      ]

However, the compiled function fails with because f[init] is held unevaluated as passed into FindRoot, I believe. Is there a way around this or another solution I am not thinking of?

You're going to have some trouble with this one no matter what, since FindRoot returns a rule as an answer and Compile only really handles functions that return numbers. Also, the slow part is almost certainly FindRoot, so Compile is unlikely to help too much. — Pillsy, Apr 05 '12 at 19:52

score 8 · Answer 1 · edited Apr 13 '17 at 12:55

I think you're going after it the wrong way. FindRoot is not compilable, and it's expected to be the most CPU-expensive part of your loop, so the possible benefits of compilation seem scarce.

To quote this most excellent answer by halirutan:

To summarize: Usually you don't need any list of supported functions, because the rule of thumb is, that compile will not work with already optimized, complicated Mathematica-methods. This includes NIntegrate, FindRoot or NMinimize. Nevertheless, Compile can easily be used to make those function-calls really fast. What you have to do is to compile your target function, because the most time with stuff like NIntegrate is spent, evaluating the integrand. The same is true for FindRoot, NMinimize and many more methods.

(emphasis mine).

score 1 · Answer 2 · edited Apr 06 '12 at 13:29

1

You should take a good look at the advanced documentation for NDSolve. NDSolve allows you to choose from many methods and you can tweak them in many ways. I think it will be very hard to beat the performance that is possible with the built in NDSolve.

edited Apr 06 '12 at 13:29

F'x

10,817
3
52
92

answered Apr 06 '12 at 00:19

Ted Ersek

7,124
19
41

Here I'm using a MWE and not the real equations. The real equations are a set of Differential-Algebraic Equations that have to be solved implicitly. I actually sent a Premier Service help request to Wolfram and they said the way I was solving it was probably the best way. I was just trying to get some performance increase here. It doesn't look like that is possible. – johntfoster Apr 06 '12 at 01:49
In that case you still might be able to improve performance if your compiled code includes a root-finding method built from scratch. It looks like you are doing root finding in 1D. If you know you can start with two points where the function changes sign Brent's method is a great choice. I could provide a Mma implementation of Brent's method. If you can't use Brent's method, another method such as Newton's method or secant method might work well. They are both easy to implement. You may also need to roll your own version of other function that can't be compiled. – Ted Ersek Apr 06 '12 at 17:38

Expression evaluation inside of FindRoot inside a Compiled Function

2 Answers2

Linked