## Re: [Vxl-users] vnl_levenberg_marquardt / gradf

 Re: [Vxl-users] vnl_levenberg_marquardt / gradf From: Mathieu Malaterre - 2007-05-18 07:50:42 ```On 5/17/07, Neil Muller wrote: > On 5/16/07, Mathieu Malaterre wrote: > > void gradf(vnl_vector const& x, vnl_matrix &J) { > > for (unsigned int i=0; i > J[0][i] = compute_a(i, x(0), x(1), x(2) ); > > } > > for (unsigned int i=0; i > J[1][i] = compute_b(i, x(0), x(1), x(2) ); > > } > > for (unsigned int i=0; i > J[2][i] = compute_c(i, x(0), x(1), x(2) ); > > } > > } > > The matrix you're calculating is not what the optimizer expects. Here, > things work much better if I replace this with: > > void gradf(vnl_vector const& x, vnl_matrix &J) { > for (unsigned int i=0; i J(i,0) = compute_a(i, x(0), x(1), x(2) ); > } > std::cout << J << std::endl; > for (unsigned int i=0; i J(i,1) = compute_b(i, x(0), x(1), x(2) ); > } > for (unsigned int i=0; i J(i,2) = compute_c(i, x(0), x(1), x(2) ); > } > } > > Personally, I find it much easier to devel with the M(i,j) syntax, > since the asserts will help you detect those sorts of errors, and only > replace those calls with M[i][j] in the final stages if it's really > performance critical. Thanks a bunch Neil ! I should have thought of transposing the matrix (actually I did not know about the '()' notation). So for later reference, here is the correct solution: void gradf(vnl_vector const& x, vnl_matrix &J) { for (unsigned int i=0; i

 [Vxl-users] vnl_levenberg_marquardt / gradf From: Mathieu Malaterre - 2007-05-16 16:09:06 ```Hi there, Alright, I think I finally understood how to use the vnl_levenberg_marquardt API to solve my problem. I am trying to fit a model function: a * ( 1. - b * exp( -t / c ) ) to a set of 8 discretized signal values (so far nothing exceptional). Here is how I formulated my problem (*). It works quite well when I am *NOT* using the gradient (gradf virtual function). But it behaves incorrectly when I start using it. Could someone please comment on my code, and let me know if there is an obvious error ? thanks for your help, -Mathieu PS: Included is a function test_deriv, I was trying to check that the value of my algebraic version of derivatives were consistant with finite difference and AFAIK it is. (*) #include #include #include #include static const double signal[] = { -13.5942, -10.2358, -7.54658, -5.39323, -3.66897, -2.28828, -1.18272, -0.297451, }; struct vnl_rosenbrock : public vnl_least_squares_function { static const unsigned int nsignals = sizeof (signal) / sizeof (double); vnl_rosenbrock(bool with_grad): vnl_least_squares_function(3, nsignals, with_grad ? use_gradient : no_gradient) {} // helpers: static double compute(double t, double a, double b, double c) { return a * ( 1. - b * exp( -t / c ) ); } static double compute_a(double t, double a, double b, double c) { return ( 1. - b * exp( -t / c) ); } static double compute_b(double t, double a, double b, double c) { return -a * exp( -t / c ); } static double compute_c(double t, double a, double b, double c) { return - t * a * b / (c * c) * exp( -t / c ); } void f(vnl_vector const& x, vnl_vector& y) { for (unsigned int i=0; i const& x, vnl_matrix &J) { for (unsigned int i=0; i myj(3,8); vnl_vector myx(3); myx[0] = myx[1] = myx[2] = valref; f.gradf(myx, myj); std::cout << "my:" << myj[0][row] << std::endl; std::cout << "my:" << myj[1][row] << std::endl; std::cout << "my:" << myj[2][row] << std::endl; vnl_vector myf_h1(8); vnl_vector myf_h2(8); double comp[3]; myx[0] = myx[1] = myx[2] = valref; myx[0] = valref - h; f.f(myx, myf_h1); myx[0] = valref + h; f.f(myx, myf_h2); comp[0] = (myf_h2[row] - myf_h1[row]) / (2 * h); myx[0] = myx[1] = myx[2] = valref; myx[1] = valref - h; f.f(myx, myf_h1); myx[1] = valref + h; f.f(myx, myf_h2); comp[1] = (myf_h2[row] - myf_h1[row]) / (2 * h); myx[0] = myx[1] = myx[2] = valref; myx[2] = valref - h; f.f(myx, myf_h1); myx[2] = valref + h; f.f(myx, myf_h2); comp[2] = (myf_h2[row] - myf_h1[row]) / (2 * h); std::cout << "comp:=" << comp[0] << std::endl; std::cout << "comp:=" << comp[1] << std::endl; std::cout << "comp:=" << comp[2] << std::endl; } int main(int, char *[]) { //test_deriv(3); vnl_rosenbrock f(true); vnl_double_3 x0(1.,1.,1.); vnl_levenberg_marquardt lm(f); vnl_vector x1 = x0.as_vector(); if (f.has_gradient()) lm.minimize_using_gradient(x1); else lm.minimize_without_gradient(x1); lm.diagnose_outcome(std::cout); std::cout << "x1 = " << x1 << std::endl; vnl_double_3 solution(3.26, 5.17, 4.5); double err = (x1 - solution).magnitude(); std::cout << "err = " << err << std::endl; return 0; } ```
 Re: [Vxl-users] vnl_levenberg_marquardt / gradf From: Neil Muller - 2007-05-17 12:51:24 ```On 5/16/07, Mathieu Malaterre wrote: > void gradf(vnl_vector const& x, vnl_matrix &J) { > for (unsigned int i=0; i J[0][i] = compute_a(i, x(0), x(1), x(2) ); > } > for (unsigned int i=0; i J[1][i] = compute_b(i, x(0), x(1), x(2) ); > } > for (unsigned int i=0; i J[2][i] = compute_c(i, x(0), x(1), x(2) ); > } > } The matrix you're calculating is not what the optimizer expects. Here, things work much better if I replace this with: void gradf(vnl_vector const& x, vnl_matrix &J) { for (unsigned int i=0; i
 Re: [Vxl-users] vnl_levenberg_marquardt / gradf From: Mathieu Malaterre - 2007-05-18 07:50:42 ```On 5/17/07, Neil Muller wrote: > On 5/16/07, Mathieu Malaterre wrote: > > void gradf(vnl_vector const& x, vnl_matrix &J) { > > for (unsigned int i=0; i > J[0][i] = compute_a(i, x(0), x(1), x(2) ); > > } > > for (unsigned int i=0; i > J[1][i] = compute_b(i, x(0), x(1), x(2) ); > > } > > for (unsigned int i=0; i > J[2][i] = compute_c(i, x(0), x(1), x(2) ); > > } > > } > > The matrix you're calculating is not what the optimizer expects. Here, > things work much better if I replace this with: > > void gradf(vnl_vector const& x, vnl_matrix &J) { > for (unsigned int i=0; i J(i,0) = compute_a(i, x(0), x(1), x(2) ); > } > std::cout << J << std::endl; > for (unsigned int i=0; i J(i,1) = compute_b(i, x(0), x(1), x(2) ); > } > for (unsigned int i=0; i J(i,2) = compute_c(i, x(0), x(1), x(2) ); > } > } > > Personally, I find it much easier to devel with the M(i,j) syntax, > since the asserts will help you detect those sorts of errors, and only > replace those calls with M[i][j] in the final stages if it's really > performance critical. Thanks a bunch Neil ! I should have thought of transposing the matrix (actually I did not know about the '()' notation). So for later reference, here is the correct solution: void gradf(vnl_vector const& x, vnl_matrix &J) { for (unsigned int i=0; i