Confusion about determining the jacobian in a rootfinding algorithm

Question

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?

Answer 1

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.