High precision multidimensional Newtons method with mpmath.findroot in Python
Question: I am trying to solve a system of equations numerically to high precision using
the multidimensional newtons method of mpmath.findroot. Here's an example
system:

    
    
    def f(x_0, x_1, x_2, x_3, x_4, x_5, y_0, y_1, y_2, y_3, y_4, y_5, l_0, l_1,l_2, l_3, l_4, l_5, l_6, l_7, l_8, l_9, l_10, l_11, l_12, l_13, l_14):
        return -x_0, -y_0 + 1, -x_5, -y_5, -x_1 - x_2, -y_1 + y_2, -x_3 - x_4, -y_3 + y_4, -(x_1 - x_5)^2 - (y_1 - y_5)^2 + 1, -(x_1 - x_4)^2 - (y_1 - y_4)^2 + 1, -(x_0 - x_5)^2 - (y_0 - y_5)^2 + 1, -(x_3 - x_4)^2 - (y_3 - y_4)^2 + 1, -(x_2 - x_3)^2 - (y_2 - y_3)^2 + 1, -(x_2 - x_5)^2 - (y_2 - y_5)^2 + 1, -(x_0 - x_5)^2 - (y_0 - y_5)^2 + 1, 2*l_10*(x_0 - x_5) + 2*l_14*(x_0 - x_5) + l_0 - .5*y_1 + .5*y_2, 2*l_9*(x_1 - x_4) + 2*l_8*(x_1 - x_5) + l_4 + .5*y_0 - .5*y_3, 2*l_12*(x_2 - x_3) + 2*l_13*(x_2 - x_5) - .5*y_0 + .5*y_4, -2*l_12*(x_2 - x_3) + 2*l_11*(x_3 - x_4) + l_6 + .5*y_1 - .5*y_5, -2*l_9*(x_1 - x_4) - 2*l_11*(x_3 - x_4) - .5*y_2 + .5*y_5, -2*l_10*(x_0 - x_5) - 2*l_14*(x_0 - x_5) - 2*l_8*(x_1 - x_5) - 2*l_13*(x_2 - x_5) + l_2 + .5*y_3 - .5*y_4, 2*l_10*(y_0 - y_5) + 2*l_14*(y_0 - y_5) + l_1 + .5*x_1 - .5*x_2, 2*l_9*(y_1 - y_4) + 2*l_8*(y_1 - y_5) + l_5 - .5*x_0 + .5*x_3, 2*l_12*(y_2 - y_3) + 2*l_13*(y_2 - y_5) + .5*x_0 - .5*x_4, -2*l_12*(y_2 - y_3) + 2*l_11*(y_3 - y_4) + l_7 - .5*x_1 + .5*x_5, -2*l_9*(y_1 - y_4) - 2*l_11*(y_3 - y_4) + .5*x_2 - .5*x_5, -2*l_10*(y_0 - y_5) - 2*l_14*(y_0 - y_5) - 2*l_8*(y_1 - y_5) - 2*l_13*(y_2 - y_5) + l_3 - .5*x_3 + .5*x_4
    import sage.libs.mpmath.all as mpmath
    startsol=[0, -0.345847373297000, 0.345847770952000, -0.499999005131000, 0.500000393924000, 0, 0.999999726908000, 0.938291180327000, 0.938290404618000, 0.404864576964000, 0.404866700310000, 0, 0, 0.0205563218420000, 0.00287356421947000, -0.0200611112418000, 0.00371301649518000, 0.00363938213640000, -0.00658839856658000, -0.00414258257348000, 0.0391290458552000, 0.162094656373000, 0.0813226022151000, 0.0974643704937000, 0.158202488927000, 0.0432804352887000, 0.0813226022151000]
    print f(*startsol)
    mpmath.mp.dps=100
    ans=mpmath.findroot(f, startsol)
    

Unfortunately it does not give me a solution but throws an error

> ZeroDivisionError: matrix is numerically singular

Is this happening because findroot is trying to calculate the jacobian? Does
this mean that the system is underdetermined?

I found the starting point using scipy's fmin_cg, but I want to polish the
solution to higher precision. As a function for fmin_cg I minimized the sum of
squares of the 27 entries of my function f.

If the problem with mpmath.findroot cannot be avoided, is there better way to
solve this system to high precision?


Answer: Yes, the system is underdetermined, its Jacobian matrix has rank 24 (at most)
instead of 27. You can check this using SymPy:

    
    
    import sympy as sp
    vars = sp.var('x_0, x_1, x_2, x_3, x_4, x_5, y_0, y_1, y_2, y_3, y_4, y_5, l_0, l_1,l_2, l_3, l_4, l_5, l_6, l_7, l_8, l_9, l_10, l_11, l_12, l_13, l_14')
    F = sp.Matrix([-x_0, -y_0 + 1, -x_5, -y_5, -x_1 - x_2, -y_1 + y_2, -x_3 - x_4, -y_3 + y_4, -(x_1 - x_5)**2 - (y_1 - y_5)**2 + 1, -(x_1 - x_4)**2 - (y_1 - y_4)**2 + 1, -(x_0 - x_5)**2 $
    J = F.jacobian(vars)
    print(J.rank())
    

Looking more carefully at the Jacobian matrix shows what the redundant
equations are. For example, these are equations numbered 0-3 and 10:

    
    
    -x_0      = 0 
    -y_0 + 1  = 0
    -x_5      = 0
    -y_5      = 0 
    -(x_0 - x_5)^2 - (y_0 - y_5)^2 + 1  = 0 
    

Clearly, the last one is a consequence of the first four. This redundancy
needs to go away.

### Using symbolic Jacobian

If mpmath.findroot is not given a function computing the Jacobian, it will use
numeric differentiation, introducing additional errors. If available, it is
better to provide a Jacobian function. Here is a small code sample for
completeness.

    
    
    import sympy as sp
    import mpmath as mp
    vars = sp.var('x, y, z')
    F = [x**2 + y**3 + z**4 - 6, x + y + z - 2, 3*x**2 + y - 5]
    J = sp.Matrix(F).jacobian(vars)
    f = lambda x0,y0,z0 : [Fc.subs(list(zip(vars, [x0,y0,z0]))) for Fc in F]
    Jac = lambda x0,y0,z0 : J.subs(list(zip(vars, [x0,y0,z0]))).tolist()
    start = [1, 1, 1]
    print(mp.findroot(f, start, J = Jac))
    

