# FindRoot: Using replacement rule inside body

I have a function similar to the following example:

f[x_, y_] := (x - 2)^2 + x ( y - 2)^2


and I would like to find the minimum with respect to x and y using FindRoot (the function is not simple enough for analytical minimization and using FindRoot on the derivatives turned out to be orders of magnitude faster than NMinimze). So far, I use

FindRoot[{D[f[x, y], x] == 0 ,
D[f[x, y], y] == 0}, {{x, 1}, {y, 3}}]
(* {x -> 2., y -> 2.} *)


which works most of the time. However, sometimes, my starting guess for x and y is off, like:

FindRoot[{D[f[x, y], x] == 0 ,
D[f[x, y], y] == 0}, {{x, 3}, {y, 20}}]
(* {x -> -4.73317*10^-30, y -> 4.} *)


yielding a wrong solution. I know that I can minimize the function for a given x reliably to obtain y as a function of x:

fy[x_] := FindRoot[{D[f[x, y], y] == 0 }, {{y, 20}}]


The question now is: Can I use this numerical function fy in the FindRoot for x?

I have tried:

FindRoot[{D[f[x, y], x] == 0 /. fy[x]}, {{x, 2}}]


which gives errors like

FindRoot::nlnum: The function value {36. x} is not a list of numbers with dimensions {1} at {y} = {20.}.


Similarly, this also does not work:

FindRoot[{D[f[x, fy[x]], x] == 0 /. fy[x]}, {{x, 2}}]
FindRoot[{D[f[x, fy[x]], x] == 0 /. fy[x]}, {{x, 2}}]


Is there a syntactically proper way to do what I intend to do?

• Can you explain why you are substituting xs for x and then back in the same equation? D[f[xs, y], xs] == 0 /. xs -> x and D[f[x, y], x] == 0 are exactly the same (D[f[xs, y], xs] == 0 /. xs -> x) === (D[f[x, y], x] == 0) gives True Commented Nov 22, 2016 at 21:53
• True, that was some extra layer of complication. I have edited the question accordingly. Commented Nov 22, 2016 at 22:32
• I still wonder though: NArgMin applied to your function returns the desired result in 23 ms (check RepeatedTiming[NArgMin[f[x, y], {x, y}]]). Is that really unacceptably slow, given that it will reliably return minima and not generic extrema as your current approach? Commented Nov 23, 2016 at 17:58
• The function provided is just an example to illustrate the issue. It is deliberately simple because I wanted to know if generally one can use FindRoot to minimize wrt. one variable first and then wrt. to the other one, without discussions on how to optimize the minimization of the specific function. In my actual application, the function is more than a page long and involves elliptical functions. FindRoot finds the minimum in a few 100 ms and NMinimize takes several minutes. Now the difference matters. Commented Nov 23, 2016 at 18:13

UPDATE: Avoiding the use of an analytical solution based on rapid brainstorming with Felix:

We have a function

f[x_, y_] := (x - 2)^2 + x (y - 2)^2


We can then take two partial derivatives

dx[x_, y_] := D[f[x, y], x]
dy[x_, y_] := D[f[x, y], y]


and we can solve for the min x at each slice of y. FindRoot is a numerical solver and it can't take symbolic arguments, so we can't get a general solution in terms of x. However, we can use the numerical solution of a partial derivative w.r.t. y == 0 to express y as a function of numeric x (when a specific x value will be plugged in by the FindRoot iterator during evaluation) and then plug it in the equation for partial derivative w.r.t. x == 0:

yi = 20;
xi = 10;
fy[x_?NumericQ] := y /. FindRoot[dy[x, y] == 0, {y, yi}]
dxnoy[x_] := dx[x, fy[x]];
{xm = x /. FindRoot[{dxnoy[x] == 0}, {{x, xi}}], fy[xm]}

(* {2., 2.} *)


Felix suggested using Block when defining partial derivatives to improve performance, see his comment below

• Yes, essentially. Just use FindRoot for xmin and I'm happy (because I can't use analytical methods in my actual application). That is, dx[x_, y_] := D[f[x, y], x] dy[x_, y_] := D[f[x, y], y] fy[x_?NumericQ] := y /. FindRoot[dy[x, y] == 0, {y, 20}] dxnoy[x_] := dx[x, fy[x]] works. If you can incorporate that in your answer I will accept it. Commented Nov 22, 2016 at 22:32
• Actually, I just realize the solution (with FindRoot for xmin) is about 100x slower than the original FindRoot[{D[f[x, y], x] == 0 , D[f[x, y], y] == 0}, {{x, 1}, {y, 3}}]. Any idea of how it comes to this performance difference? Commented Nov 22, 2016 at 22:54
• @Felix I see, this is probably not going to work in general case -- FindRoot is a numerical solver, so symbol x tells nothing and it complains. In your specific case of this function, you can see that dy[x, y] is 2 x (-2 + y), so it only is equal to 0 when either y=2 or x=0. This gives you a solution immediately, and you don't even need to do a second FindRoot. Commented Nov 22, 2016 at 22:57
• Well, the function f is just an example. My actual function is a page long and confidential. What matters is that f is analytical and so are it's derivatives, but finding the roots of the derivatives can't be done analytically. Commented Nov 22, 2016 at 23:01
• I found the performance issue with my previous suggestion based on your code. This works with excellent timing: dx[x_, y_] := Block[{xs}, D[f[xs, y], xs] /. xs -> x]; dy[x_, y_] := Block[{ys}, D[f[x, ys], ys] /. ys -> y]; fy[x_?NumericQ] := y /. FindRoot[dy[x, y] == 0, {y, 20}]; dxnoy[x_] := dx[x, fy[x]];FindRoot[{dxnoy[x] == 0}, {{x, 5}}]. It is important to take the derivative wrt. a variable, not a number. Hence, calling dx[x_, y_] := D[f[x, y], x] with some number for x doesn't really work. With dx[x_, y_] := Block[{xs}, D[f[xs, y], xs] /. xs -> x] the problem does not exist. Commented Nov 22, 2016 at 23:12

For this you can use FindRoots2D. With your f[x,y] and the intervall [-10,10]:

roots = FindRoots2D[Grad[f[x, y], {x, y}], {x, -10, 10}, {y, -10, 10}] // Chop
{{0, 0}, {0, 4.}, {2., 2.}}


Now we search for the minimal point:

{{x, y}, f} = Flatten[MinimalBy[Thread@{roots, f @@@ roots}, Last], 1]
{{2., 2.}, 0.}


FindRoots2D you can find here

I test the solution analytical:

Clear[x, y, f]
f[x_, y_] = (x - 2)^2 + x (y - 2)^2;
pts = Solve[Grad[f[x, y], {x, y}] == 0, {x, y}]
{{x -> 0, y -> 0}, {x -> 2, y -> 2}, {x -> 0, y -> 4}}

hesse = D[f[x, y], {{x, y}, 2}] /. pts;
{PositiveDefiniteMatrixQ /@ hesse, NegativeDefiniteMatrixQ /@ hesse} // Column
{False, True, False}, {False, False, False}


Point {2,2} is the minimum

• I think this is fundamentally the same as my approach as mine, isn't it? I mean in the FindRoots2D code it reads FindRoot[funcs, {x, #}, {y, #2}, seq]. Of course, all this ContourPlot stuff is very fancy but it is too slow for my application. For my present example, I would really like to use the fact that f can be minimized wrt. y for any given x. Commented Nov 22, 2016 at 21:28
• @Felix So far you didn't make a mistake. You have only found 2 extremal points {2,2} and {0,4} by try and you are terrified for this result. But you have not found all extremal points. For this you can use FindRoot2D with advantage. If you study this extremly valuable procedure, you see how to find the seeds (starting points) for FindRoot for which you are looking for. it makes no sense to follow a path wich does not lead to the goal!
– user36273
Commented Nov 22, 2016 at 22:01
• I don't say your solution is wrong. It is just too time consuming. I have to find minima for millions of such functions and performance really matters. In the present example, the originally proposed FindRoot[{D[f[x, y], x] == 0 , D[f[x, y], y] == 0}, {{x, 1}, {y, 3}}] finds the solution in less than 1 ms. I am aiming for something of similar speed, but more robust. Commented Nov 22, 2016 at 22:59
• @Felix But this is not a minima! If you have found a systematic way to find all roots with FindRoot in a shorter time than FindRoots2D, without knowing how many and where they are located, please let it me know. Try it with a just slightly modified function f[x_,y_]=(x-2)^2+x (y-2)^2+Cos[x y]^2. On my notebook I need 0.39 s.
– user36273
Commented Nov 23, 2016 at 10:39
• I agree, for finding all roots your suggestion is probably best. I should have mentioned explicitly: I know that the function has at most one minimum. Also, 0.39s should be compared to ~1ms for the accepted solution. Commented Nov 23, 2016 at 15:27