# Confusion about determining the jacobian in a rootfinding algorithm

I have written some Python code to determine the numerical roots of the following non-linear equation: $$f_m=\tan\lambda_m - \frac{\lambda_m}{1+a}$$ where $\lambda_m\gt0$ and $a\geq0$. The code is:

import numpy as np
from scipy.optimize import root

func = lambda λm, a: np.tan(λm) - λm/(1+a)
λm0 = np.linspace(np.pi, m*np.pi, m) # initial guess
sol = root(func, λm0, args=(a), method='lm') # least-square optimization


This works fine, e.g. for input $m, a = 3, 1$ the result is: $$\lambda_m=[4.27478227, 7.59654602, 10.81267333]$$ which are the roots to $f_m=0$ for $m=1,2,3$ as can be checked easily.

As a next step, I would like to provide the Jacobian, which i think should be: $$J_{m,n}=\frac{\partial f_{m}}{\partial\lambda_{n}}=\begin{cases} \frac{1}{\cos^2\lambda_{m}}-\frac{1}{1+a} & m=n\\ 0 & m\neq n \end{cases}$$ i.e. in matrix-form with the above values for $\lambda_m$: $$\boldsymbol{J}=\begin{bmatrix} 5.06844087 & 0 & 0\\ 0 & 14.92687786 & 0\\ 0 & 0 & 29.72847616 \end{bmatrix}$$

However, the approximated Jacobian being output by the optimization routine is: $$\boldsymbol{\tilde{J}}=\begin{bmatrix} -29.72850523 & 0 & 1\\ 0 & -14.92688569 & 0\\ 0 & 0 & 5.06844398 \end{bmatrix}$$

Clearly the magnitudes of approximated derivatives on the diagonal correspond to my analytical values for the Jacobian, but the order of the values is somehow different (running backwards from $n=3,2,1$) and some of the values are negative. Furthermore, an additional $\tilde{J}_{1,3}=1$ is present (more of these spurious ones seem to appear as $n$ grows) which shouldn't be present.

Can someone explain the differences between the approximate and analytical Jacobian?

• Sorry, stupid question. In your non-linear equation, what is the independent variable you want to solve for? and what the parameters? what is $f_m$? and what is $m$? – Walter Jun 7 '17 at 18:54
• @Walter the independent variable is $\lambda_m$, $\alpha$ is just a constant, $f_m$ is the equation of which I want to find the roots. $m$ takes integer values and represents the fact that there are infinite solutions to the non-linear equation, i.e. $\lambda_m$ is an infinitely long vector. Of course I can't represent this numerically so I choose a finite number of solutions, e.g. $n=3$ in the above example, but it will be scaled up to several thousands of solutions. – nluigi Jun 7 '17 at 21:29
• @nluigi You know this is a one-dimensional problem, right? (Just to check—I'm assuming you're asking about the Jacobian and you aren't actually solving a 1d equation like this.) – Kirill Jun 7 '17 at 21:59
• From scipy.optimize I tried brentq (like this: brentq(f, (k - 0.5)*pi + 1e-8, (k + 0.5)*pi - 1e-8)), and it computed the first 1000 roots in 25ms. This seems to work up to $k\leq 10^7$, but you could probably come up with a better bracket with a little more work. – Kirill Jun 8 '17 at 18:59
• I think if you write out what lmdif is doing (Levenberg-Marquardt), you'll find that it's applying basically Newton's method to your 1-d problem, in parallel to all the points at once, but with the extra work that comes from not knowing that the Jacobian is diagonal. And you also miss out on being able to bracket the roots. – Kirill Jun 8 '17 at 19:05

If you look at the full output of your script, sol has fields ipvt and qtf. Both of these fields come from MINPACK (http://www.netlib.org/minpack/lmdif.f), which is the backend here (https://github.com/scipy/scipy/blob/2526df72e5d4ca8bad6e2f4b3cbdfbc33e805865/scipy/optimize/_root.py#L246; https://docs.scipy.org/doc/scipy-0.19.0/reference/generated/scipy.optimize.leastsq.html).

The first, ipvt explains what happened to variable ordering—there was a pivoting step when dealing with the Jacobian and when I ran your code, the pivots were $[3,2,1]$, matching the permutation of the matrix. The second, qft is given in the documentation as (q transpose)*fvec, which in my case is [ -2.52376384e-07, -1.46294147e-07, 5.03946448e-07]. This explains the sign changes for the first two entries, but the magnitudes are off by $10^7$—I don't know why.

One other thing I noticed is that lmdif doesn't produce the actual Jacobian—note that the field is called fjac (as in lmdif), not jac (as in scipy). Quote:

c       fjac is an output m by n array. the upper n by n submatrix
c         of fjac contains an upper triangular matrix r with
c         diagonal elements of nonincreasing magnitude such that
c
c                t     t           t
c               p *(jac *jac)*p = r *r,
c
c         where p is a permutation matrix and jac is the final
c         calculated jacobian. column j of p is column ipvt(j)
c         (see below) of the identity matrix. the lower trapezoidal
c         part of fjac contains information generated during
c         the computation of r.


I haven't tested this, but there would probably be even more inconsistency if you tried a coupled system of equations.

• Ok, so basically fjac is not the jacobian... that explains the difference. Thx! – nluigi Jun 8 '17 at 7:53