implement findroot(lambda x, y: x+y, (1, 1))

[GoogleCodeExporter]

In r767 I implemented a multidimensional nonlinear solver. What do you
think how findroot() should be called in the multidimensional case?

Assuming something like

def f(x1, ..., xn):
    return (f1(x1, ..., xn), ..., fm(x1, ..., xn))

How should findroot() recognize it's a multidimensional function?

A possibility is evaluating f(*x0) an see whether the result is scalar. Or
there could be a flag. Or f is a list of functions and findroot()
transforms it to a function like defined above.

What do you think?

Original issue reported on code.google.com by Vinzent.Steinberg@gmail.com on 19 Nov 2008 at 11:32

[GoogleCodeExporter]

Evaluating f(*x0) seems robust enough to me.

It would be convenient if you could pass a list of functions, so I think this 
should
be supported. But it is not sufficient, as the components can not generally be
computed independently.

Original comment by fredrik....@gmail.com on 19 Nov 2008 at 11:40

[GoogleCodeExporter]

Evaluating this would mean an unnecessary evaluation (if it's not reused which 
would
add much complexity). But there is no other way, except for letting the user 
specify
it explicitly. And the performance impact is marginal. So let's do it this way.

Original comment by Vinzent.Steinberg@gmail.com on 19 Nov 2008 at 12:07

[GoogleCodeExporter]

r722

Original comment by Vinzent.Steinberg@gmail.com on 19 Nov 2008 at 9:06

• Changed state: Fixed

[GoogleCodeExporter]

This should work: findroot(lambda x, y: x+y, (1, 1))

Original comment by fredrik....@gmail.com on 19 Nov 2008 at 9:14

[GoogleCodeExporter]

The problem is that the method to evaluate f(*x0) does not work here. It should
however work with multidimensional=True, which it does not. Is there a way to 
get the
number of arguments of f?

Original comment by Vinzent.Steinberg@gmail.com on 21 Nov 2008 at 2:55

• Changed title: _fix findroot(lambda x, y: x+y, (1, 1)) _
• Changed state: Started

[GoogleCodeExporter]

>>> f = lambda x, y: x+y
>>> f(*(1,1))
2

What is the problem?

Original comment by fredrik....@gmail.com on 21 Nov 2008 at 2:58

[GoogleCodeExporter]

I thinks that it's a one dimensional function and thus tries to use the secant 
method.

Original comment by Vinzent.Steinberg@gmail.com on 21 Nov 2008 at 3:01

[GoogleCodeExporter]

It should be able to deduce the dimension from x0 having length, no?

Original comment by fredrik....@gmail.com on 21 Nov 2008 at 3:09

[GoogleCodeExporter]

For the secant method x0 can be twodimensional, for other onedimesional solvers 
even
threedimensional. So this works only for len(x0) > 3.

Original comment by Vinzent.Steinberg@gmail.com on 21 Nov 2008 at 3:17

[GoogleCodeExporter]

Right.

You could use f.func_code.co_argcount, but this is not always reliable (e.g. for
functions using *args)

Original comment by fredrik....@gmail.com on 21 Nov 2008 at 3:23

[GoogleCodeExporter]

Function using *args are usually multidimensional, because you don't want to use
*args for one argument. So it should be implemented this way.

BTW, findroot(lambda x, y: x+y, (1, 1)) attempts to solve an underdetermined 
system.
Those are not solvable using Newton's method.

Original comment by Vinzent.Steinberg@gmail.com on 21 Nov 2008 at 3:35

[GoogleCodeExporter]

Doesn't being underdetermined simply mean that you can pick an arbitrary 
solution?

Original comment by fredrik....@gmail.com on 21 Nov 2008 at 3:56

[GoogleCodeExporter]

This would also be a possibility:

findroot(f, (x, y)) -- solve 2D problem with start vector (x,y)
findroot(f, x1, x2) -- solve problem with two start vectors x1, x2 given to 
solver

Original comment by fredrik....@gmail.com on 21 Nov 2008 at 3:59

[GoogleCodeExporter]

Being underdetermined means that the solution is a function. You can transform 
x+y==0
to x==-y and you got a solution depending on another unknown. However you have 
turn
some variables into constants.

Or do you mean that one solution should be returned?

However, the current solver does not support this, so this is not an issue. 
Currently
findroot(lambda x, y: x+y, (1, 1)) should not work. Or how would you implement 
it?

Original comment by Vinzent.Steinberg@gmail.com on 22 Nov 2008 at 12:16

• Changed title: _implement findroot(lambda x, y: x+y, (1, 1)) _
• Changed state: New

[GoogleCodeExporter]

> Or do you mean that one solution should be returned?

Yes, i think it would be useful to return an arbitrary solution.

Original comment by fredrik....@gmail.com on 22 Nov 2008 at 9:47

[GoogleCodeExporter]

On a related note:

>>> findroot(lambda x: 0, 1)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "c:\source\mp\trunk\mpmath\settings.py", line 135, in g
    return f(*args, **kwargs)
  File "c:\source\mp\trunk\mpmath\optimization.py", line 835, in findroot
    if not isinstance(x, (list, tuple, matrix)):
UnboundLocalError: local variable 'x' referenced before assignment

Original comment by fredrik....@gmail.com on 22 Nov 2008 at 9:50

[GoogleCodeExporter]

> Yes, i think it would be useful to return an arbitrary solution.

So let's assign random values to random variables until we get a determined 
system?
Shouldn't this rather be done by the user?

I think findroot(lambda x: 0, 1) fails because the generator does not return
anything. Probably it cancels immediately due to zero division, assuming the 
previous
approximation has already been yielded, which is no true in this case.

Original comment by Vinzent.Steinberg@gmail.com on 22 Nov 2008 at 1:04

[GoogleCodeExporter]

I think findroot should try to find *a* root, even it's trivial or not uniquely
defined.  (Even in the "determined" case, the root is not generally uniquely
determined globally.) And if f(*x0) is zero, x0 should be returned with no 
further
questions asked :)

How do I as a user know which random values of x and y give a root of f(x,y)? 
It's
easy if I know that f(x,y) = x+y, but in general f might be something 
complicated
that just happens to be underdetermined by chance. If it's easy for the 
rootfinder to
pick a zero, then it should do it.

Original comment by fredrik....@gmail.com on 22 Nov 2008 at 2:17

[GoogleCodeExporter]

I fixed findroot() for f(x0) == 0.

> If it's easy for the rootfinder to
> pick a zero, then it should do it.

How? Let's say we have f(a,b,c,d) == e. Choose a=b=c=0 and solve for d? I don't 
want
the solver to do this, IMHO the user should care about such things.

For lambda x,y: x+y it would make more sense to compute the Hessian matrix to
minimize the function.

Original comment by Vinzent.Steinberg@gmail.com on 23 Nov 2008 at 5:29

[GoogleCodeExporter]

You can just use Newton's method as usual, no?

Only each time you solve the Jacobian equation, there are infinitely many 
solutions
and you can pick one arbitrarily.

Consider my example f(x,y) = x+y+1 with (x0,y0) = (1,1). At the first step you 
get
the equation [1,1]*[x1-1,y1-1]^T = -f(x0,y0) = -3 with infinitely many 
solutions x1 =
-1-y1. Pick, say, y1 = 0. This gives (-1, 0) which turns out to be an exact 
root of f
in this case (because f is linear). In general, with f being nonlinear, you've 
ended
up closer to the root. Repeat.

> How? Let's say we have f(a,b,c,d) == e. Choose a=b=c=0 and solve for d?

This does not work because a function might not have zeros for all fixed 
choices.
Consider f(x,y) = (x-1)^2 + y^2. This has solution for x = 1 but not for any 
other
choice of x (assuming real numbers).

An alternative solution is to add a dummy output variable, e.g. either return
(f(x,y), f(x,y)) or a constant zero (f(x,y), 0). Mathematica can solve (x-1)^2 
+ y^2
if you do this:

In[7]:= FindRoot[(x-1)^2+y^2, {{x, 2}, {y, 2}}]

FindRoot::nveq:
   The number of equations does not match the number of variables in
                    2    2
    FindRoot[(x - 1)  + y , {{x, 2}, {y, 2}}].

                        2    2
Out[7]= FindRoot[(x - 1)  + y , {{x, 2}, {y, 2}}]

In[8]:= FindRoot[{(x-1)^2+y^2, 0}, {{x, 2}, {y, 2}}]

                                 -9
Out[8]= {x -> 1., y -> 7.45058 10  }

But I don't see any point in requiring workarounds of this kind. Either way, the
solver ends up with a Jacobian equation with less than full rank, and has to 
know how
to deal with that.

I do think findroot "just work", in the most natural way, for an input like 
(x-1)^2 +
y^2.

Original comment by fredrik....@gmail.com on 23 Nov 2008 at 6:00

[GoogleCodeExporter]

Well, this should work:

>>> findroot([lambda x,y:(x-1)**2+y**2, lambda x,y:0], (2, 2))
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "mpmath/settings.py", line 135, in g
    return f(*args, **kwargs)
  File "mpmath/optimization.py", line 827, in findroot
    for x, error in iterations:
  File "mpmath/optimization.py", line 561, in __iter__
    s = lu_solve(Jx, fxn)
  File "mpmath/settings.py", line 135, in g
    return f(*args, **kwargs)
  File "mpmath/linalg.py", line 209, in lu_solve
    A, p = LU_decomp(A)
  File "mpmath/linalg.py", line 136, in LU_decomp
    * absmin(A[k,j])
  File "mpmath/mptypes.py", line 224, in __rdiv__
    return make_mpf(mpf_rdiv_int(t, s._mpf_, prec, rounding))
  File "mpmath/libmpf.py", line 822, in mpf_rdiv_int
    return mpf_div(from_int(n), t, prec, rnd)
  File "mpmath/libmpf.py", line 781, in mpf_div
    raise ZeroDivisionError
ZeroDivisionError
>>> findroot([lambda x,y:(x-1)**2+y**2]*2, (2, 2))
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "mpmath/settings.py", line 135, in g
    return f(*args, **kwargs)
  File "mpmath/optimization.py", line 827, in findroot
    for x, error in iterations:
  File "mpmath/optimization.py", line 561, in __iter__
    s = lu_solve(Jx, fxn)
  File "mpmath/settings.py", line 135, in g
    return f(*args, **kwargs)
  File "mpmath/linalg.py", line 211, in lu_solve
    x = U_solve(A, b)
  File "mpmath/linalg.py", line 182, in U_solve
    x[i] /= U[i,i]
  File "<string>", line 8, in __div__
  File "mpmath/libmpf.py", line 777, in mpf_div
    if t == fzero: raise ZeroDivisionError
ZeroDivisionError

Another way to fix it is to support underdetermined systems in lu_solve.

Original comment by Vinzent.Steinberg@gmail.com on 2 Dec 2008 at 9:46

[GoogleCodeExporter]

Yes, the problem seems to be that lu_solve divides by zero when a zero row 
turns up.
With this fixed, one could just handle the underdetermined case by appending 
dummy
zero rows.

I just found another instance where this is causing trouble: if you apply 
pade() to a
rational function, it fails if the denominator polynomial has higher degree 
than the
actual function. For example with f(x) = 1/(1-x), pade([1,1,1,1,1], 3,1) works 
but
pade([1,1,1,1,1,1], 3,2) doesn't.

So let's fix this (I can do it eventually, but I have other issues to deal with 
right
now).

Original comment by fredrik....@gmail.com on 3 Dec 2008 at 12:08