I've got an issue when I use vnl_levenberg_marquardt to refine the estimation of the 2-D homography between calibration rig plane and image plane.
My algorithm follows the gold standard algorithm proposed in [Hartley & Zisserman, 2000] (algorithm 3.3, p98) and I choose to use Sampson error for the sake of speed. For implementation, I write a class similar to FMatrixComputeNonLinear. To gain faster convergence, I provide a customized routine gradf(), in which the Jocabian is a 2N x 9 matrix (equation 3.31, p114), where N is the number of point correspondences between the two planes. 
Here comes the problem. For geometric minimization of Sampson error, the number of residuals is N, since there is one residual (the norm $||\delta_x||^2$, equation 3.12, p83) for each correspondence:
HMatrix2DComputeSampson::HMatrix2DComputeSampson(int N)
: vnl_least_square_function(9, N, use_gradient) {}
However, if my class is initialized like this, I find in the vnl_levenberg_marquardt code that dimension of Jocabian matrix in gradf() is fixed as N x 9 and my program will crash. One workaround is to set the number of residuals to 2*N.
HMatrix2DComputeSampson::HMatrix2DComputeSampson(int N)
: vnl_least_square_function(9, 2*N, use_gradient) {}
But this will introduce some minor problems (more RAM consumption, slower speed, unintuitive).
In all, I'd like to know your opinions on this issue.
I'd appreciate your patience to read this mail and thanks in advance.
Lianqing Yu
PhD Candidate
Robot Vision Group
National Laboratory of Pattern Recognition (NLPR)
Institute of Automation, Chinese Academy of Sciences (CASIA)
No. 95, Zhongguancun East Road,
100080,Beijing, China
Email: lqyu@nlpr.ia.ac.cn