Making the Jacobian analytical brings some gain. However, Newton does not converge when the theta is larger than around 1.25 for the parameters in #2802. The loss landscape is very bad for Newton convergence, that leads to extreme oscillations.
With the previous IterativeUndistortion the oscillation can lead to a different root potentially with norm larger than pi/2, which essentially invalids the tangent function, leading to > 500 pixel error. Adding a trust region would reduce the forward-backward error to always less than 75 pixels (extreme case happens for thetas close to pi/2), but does not resolve the convergence issue. The fact that k3, k4 are large enough for the radial distortion with an exponentially growing theta^8 makes it hard to converge with theta that is larger than 1.25.
I dont think we can do better with Newton here. I also tried adding damping factors that does not give much gain. I believe we need to resort to other optimization approach in order to fully resolve the issue.