3

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?

Felix
  • 3,871
  • 13
  • 33
  • 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 – Stitch Nov 22 '16 at 21:53
  • True, that was some extra layer of complication. I have edited the question accordingly. – Felix Nov 22 '16 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? – MarcoB Nov 23 '16 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. – Felix Nov 23 '16 at 18:13

2 Answers2

2

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

Community
  • 1
  • 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. – Felix Nov 22 '16 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! –  Nov 22 '16 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. – Felix Nov 22 '16 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. –  Nov 23 '16 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. – Felix Nov 23 '16 at 15:27
  • @Felix I think there's a misunderstanding how to find a minimum of a function. One have to find first all roots of the gradient of the function and then investigate using this roots on what positions the function is minimal. Finding only one root with FindRoot is not determining the minimum, you may only get a maximal, minimal or saddle point. Try my function in the comment, there are 7 roots and only one is the minimum in the intervall [-5,5], found in 0.16 s. –  Nov 23 '16 at 19:39
2

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

Stitch
  • 4,205
  • 1
  • 12
  • 28
  • 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. – Felix Nov 22 '16 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? – Felix Nov 22 '16 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. – Stitch Nov 22 '16 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. – Felix Nov 22 '16 at 23:01
  • 1
    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. – Felix Nov 22 '16 at 23:12
  • @Felix in this case, you need to do a bit more analysis of your function -- the problem is that you have a minimum {2,2} and saddle points ({0,0} and {0,4}) where both partial derivatives are 0. Your original algorithm finds those values depending on where you start. So, if you can't use FindMinimum[f[x, y], {{x, 10}, {y, 3}}] you have to consider using other numerical methods. – Stitch Nov 22 '16 at 23:17
  • Your solution answers my question, just replace Solve by FindRoot as I have suggested in my last comment. If you do this edit I will accept the answer. – Felix Nov 22 '16 at 23:27
  • @Felix great ideas, updated with reference to these comments – Stitch Nov 23 '16 at 00:08