# Custom jacobian (dfdp) in nonlin_curvefit - calling parameters?

12 messages
Open this post in threaded view
|

## Custom jacobian (dfdp) in nonlin_curvefit - calling parameters?

 Hi, I am trying to find out, if providing custom jacobian results in faster execution than the default real finite differences. However, I got stuck at defining the custom dfdp function. Please is there any example which involves multiple independent values? My case is fitting sum of two sine functions at vector of times t to observed values y. The problem is I do not know how the custom dfdp function should look like, with what params it is being called. In my code below my dfdp (f_dfdp) gets called with row of times (t) which I cannot convert to the jacobian matrix. The dfdp description in nonlin_residmin talks about some structures, but I was not able to find any example. As of the unit test in nonlin_curfit.m (chittkowski_327) - I was able to find the commit diff but this case again looks like a single independent value is being used. I very much appreciate any help, thanks a lot. Best regards, Pavel. t = 0:100; function y = f_calc(f1, a1, s1, f2, a2, s2, t)    persistent PI2 = 2 * pi;    y = a1 * cos(PI2 * f1 * t + s1) + a2 * cos(PI2 * f2 * t + s2); endfunction function J = f_dfdp(f1, a1, s1, f2, a2, s2, t)    persistent PI2 = 2 * pi;      J = [-PI2 * a1 * t * sin(PI2 * f1  * t + s1) , cos(PI2 * f1 * t + s1) ,  -a1  * sin(PI2 * f1 * t + s1), -PI2 * a2 * t * sin(PI2 * f2 * t + s2),  cos(PI2 * f2 * t + s2), -a2 * sin(PI2 * f2 * t + s2)]; endfunction f = @(p, t) f_calc(p(1), p(2), p(3), p(4), p(5), p(6), t); my_dfdp= @(p, t) f_dfdp(p(1), p(2), p(3), p(4), p(5), p(6), t); settings = optimset ('dfdp', my_dfdp, 'Display','iter'); init = [measF1; measA1; measPhase1;measF2; measA2; measPhase2]; [p, model_values, cvg, outp] = nonlin_curvefit(f, init, t, y, settings);
Open this post in threaded view
|

## Re: Custom jacobian (dfdp) in nonlin_curvefit - calling parameters?

 On Fri, Nov 15, 2019 at 02:54:22PM +0100, Pavel Hofman wrote: > Hi, > > I am trying to find out, if providing custom jacobian results in faster > execution than the default real finite differences. However, I got stuck at > defining the custom dfdp function. Please is there any example which > involves multiple independent values? > > My case is fitting sum of two sine functions at vector of times t to > observed values y. The problem is I do not know how the custom dfdp function > should look like, with what params it is being called. In my code below my > dfdp (f_dfdp) gets called with row of times (t) which I cannot convert to > the jacobian matrix. The dfdp description in nonlin_residmin talks about > some structures, but I was not able to find any example. The interface of your functions and the way you call nonlin_curvefit seem ok. The code inside f_dfdp() doesn't seem ok. I don't think the documentation of the optim package contains an example for explicitly computing a Jacobian with m-code (because this is not specific to the optim package). But in the code of 'optim_problems.m' there is an example:   ret.curve.schittkowski_327.dfdp = ...       @ (x, p) [1 + exp(-p(2) * (x - 8)), ...                 (p(1) + .49) * (8 - x) .* exp(-p(2) * (x - 8))]; Note that 'x' (or 't', in your code) is a column vector. Olaf > As of the unit test in nonlin_curfit.m (chittkowski_327) - I was able to > find the commit diff but this case again looks like a single independent > value is being used. > > I very much appreciate any help, thanks a lot. > > Best regards, > > Pavel. > > > t = 0:100; > > function y = f_calc(f1, a1, s1, f2, a2, s2, t) >   persistent PI2 = 2 * pi; > >   y = a1 * cos(PI2 * f1 * t + s1) + a2 * cos(PI2 * f2 * t + s2); > endfunction > > function J = f_dfdp(f1, a1, s1, f2, a2, s2, t) >   persistent PI2 = 2 * pi; > >     J = [-PI2 * a1 * t * sin(PI2 * f1  * t + s1) , cos(PI2 * f1 * t + s1) , > -a1  * sin(PI2 * f1 * t + s1), -PI2 * a2 * t * sin(PI2 * f2 * t + s2), > cos(PI2 * f2 * t + s2), -a2 * sin(PI2 * f2 * t + s2)]; > endfunction > > > f = @(p, t) f_calc(p(1), p(2), p(3), p(4), p(5), p(6), t); > > my_dfdp= @(p, t) f_dfdp(p(1), p(2), p(3), p(4), p(5), p(6), t); > > settings = optimset ('dfdp', my_dfdp, 'Display','iter'); > > > init = [measF1; measA1; measPhase1;measF2; measA2; measPhase2]; > > > [p, model_values, cvg, outp] = nonlin_curvefit(f, init, t, y, settings); > > > > > -- public key id EAFE0591, e.g. on x-hkp://pool.sks-keyservers.net signature.asc (849 bytes) Download Attachment
Open this post in threaded view
|

## Re: Custom jacobian (dfdp) in nonlin_curvefit - calling parameters?

 Hi Olaf, > > The interface of your functions and the way you call nonlin_curvefit > seem ok. The code inside f_dfdp() doesn't seem ok. I don't think the > documentation of the optim package contains an example for explicitly > computing a Jacobian with m-code (because this is not specific to the > optim package). But in the code of 'optim_problems.m' there is an > example: > >    ret.curve.schittkowski_327.dfdp = ... >        @ (x, p) [1 + exp(-p(2) * (x - 8)), ... > (p(1) + .49) * (8 - x) .* exp(-p(2) * (x - 8))]; > > Note that 'x' (or 't', in your code) is a column vector. Thanks a lot for your help. Function APIs of that example are: ret.curve.schittkowski_327.f = @(x, p) ...... ret.curve.schittkowski_327.dfdp = @(x, p) ...... while nonlin_curvefit has @(p, x). Should the partial derivative dfdp function for nonlin_curvefit be @(x, p), or @(p, x)? One more related question. In my other code I nonlin_curvefit 4 parameters to two equations (each using all of the params), again for a row/vector of independent times t (i.e. x). As a result, my f(p, x) returns a 2D matrix. Works perfect. When providing a custom dfdp function, what dimensions should its output have? 2D + 1D for each p, i.e. 3D? Thanks a lot, Pavel.
Open this post in threaded view
|

## Re: Custom jacobian (dfdp) in nonlin_curvefit - calling parameters?

 Dne 16. 11. 19 v 9:30 Pavel Hofman napsal(a): > Hi Olaf, > >> >> The interface of your functions and the way you call nonlin_curvefit >> seem ok. The code inside f_dfdp() doesn't seem ok. I don't think the >> documentation of the optim package contains an example for explicitly >> computing a Jacobian with m-code (because this is not specific to the >> optim package). But in the code of 'optim_problems.m' there is an >> example: >> >>    ret.curve.schittkowski_327.dfdp = ... >>        @ (x, p) [1 + exp(-p(2) * (x - 8)), ... >>         (p(1) + .49) * (8 - x) .* exp(-p(2) * (x - 8))]; >> >> Note that 'x' (or 't', in your code) is a column vector. > > Thanks a lot for your help. > > Function APIs of that example are: > > ret.curve.schittkowski_327.f = @(x, p) ...... > ret.curve.schittkowski_327.dfdp = @(x, p) ...... > > while nonlin_curvefit has @(p, x). Should the partial derivative dfdp > function for nonlin_curvefit be @(x, p), or @(p, x)? Trial/error shows the @(p, x) dfdp API is the case for nonlin_curvefit. Works perfect now when I fixed my dfdp function to use element-by-element multiplication. As of performance - The precisely calculated dfdp function never gives worse results than real finite differences. Time performance in my case depends on length of the independent vector - from 60% to 130% of the default real finite differences. I just wonder how the 2D case should be handled. Thanks a lot, Pavel.
Open this post in threaded view
|

## Re: Custom jacobian (dfdp) in nonlin_curvefit - calling parameters?

 On Sat, Nov 16, 2019 at 08:24:26PM +0100, Pavel Hofman wrote: > I just wonder how the 2D case should be handled. The Jacobian function is called with the 'x' reshaped to a column vector (like 'x(:)' would do it). The output should be calculated for this column vector (one row for each element of the column vector and one column for each partial derivative). Olaf -- public key id EAFE0591, e.g. on x-hkp://pool.sks-keyservers.net signature.asc (849 bytes) Download Attachment
Open this post in threaded view
|

## Re: Custom jacobian (dfdp) in nonlin_curvefit - calling parameters?

 Hi Olaf, Dne 17. 11. 19 v 7:52 Olaf Till napsal(a): > On Sat, Nov 16, 2019 at 08:24:26PM +0100, Pavel Hofman wrote: >> I just wonder how the 2D case should be handled. > > The Jacobian function is called with the 'x' reshaped to a column > vector (like 'x(:)' would do it). The output should be calculated for > this column vector (one row for each element of the column vector and > one column for each partial derivative). Now I have for f (using the default real finite differences): f = @(p, x) ... where x is one row of indeps, and f returns two rows, each row for one equation (each equation uses all p elements). This code works OK. If I wanted to use explicit jacobian to improve precision and calculation stability, please how should I handle the two equations in f and dfdp? Something like this? The original row of x would be transformed to column twice the length, second half duplicate of the first one. f output = column vector where first half are outputs of the first equation (for the first duplicate of the original x), and the second half are outputs of the second equation (the the second duplicate of the original x). Otherwise I do not know how to convert outputs of two equations to a single column vector of the same length as the indeps vector. I very much appreciate your help. Best regards, Pavel.
Open this post in threaded view
|

## Re: Custom jacobian (dfdp) in nonlin_curvefit - calling parameters?

 On Sun, Nov 17, 2019 at 10:32:19AM +0100, Pavel Hofman wrote: > Hi Olaf, > > Dne 17. 11. 19 v 7:52 Olaf Till napsal(a): > > On Sat, Nov 16, 2019 at 08:24:26PM +0100, Pavel Hofman wrote: > > > I just wonder how the 2D case should be handled. > > > > The Jacobian function is called with the 'x' reshaped to a column > > vector (like 'x(:)' would do it). The output should be calculated for > > this column vector (one row for each element of the column vector and > > one column for each partial derivative). > > > Now I have for f (using the default real finite differences): > > f = @(p, x) ... where x is one row of indeps, and f returns two rows, each > row for one equation (each equation uses all p elements). This code works > OK. > > If I wanted to use explicit jacobian to improve precision and calculation > stability, please how should I handle the two equations in f and dfdp? Sorry, please forget what I said about the reshaping of 'x', I was confused... I probably shouldn't try to answer help requests while I'm working on something different. And I probably don't understand what you mean with 'the 2D case'. The Jacobian should just compute one row of partial derivatives for each element of the result of 'f'. HTH, Olaf -- public key id EAFE0591, e.g. on x-hkp://pool.sks-keyservers.net signature.asc (849 bytes) Download Attachment
Open this post in threaded view
|

## Re: Custom jacobian (dfdp) in nonlin_curvefit - calling parameters?

 Hi Olaf, >> f = @(p, x) ... where x is one row of indeps, and f returns two rows, each >> row for one equation (each equation uses all p elements). This code works >> OK. >> >> If I wanted to use explicit jacobian to improve precision and calculation >> stability, please how should I handle the two equations in f and dfdp? > > Sorry, please forget what I said about the reshaping of 'x', I was > confused... I probably shouldn't try to answer help requests while I'm > working on something different. And I probably don't understand what > you mean with 'the 2D case'. > > The Jacobian should just compute one row of partial derivatives for > each element of the result of 'f'. Thanks but my f returns a 2D matrix (a row for each equation, a column for each element of the indeps vector x). Should the dfdp return a 3D matrix? Thanks a lot. Pavel.
Open this post in threaded view
|

## Re: Custom jacobian (dfdp) in nonlin_curvefit - calling parameters?

 On Mon, Nov 18, 2019 at 03:55:41PM +0100, Pavel Hofman wrote: > Hi Olaf, > > > > f = @(p, x) ... where x is one row of indeps, and f returns two rows, each > > > row for one equation (each equation uses all p elements). This code works > > > OK. > > > > > > If I wanted to use explicit jacobian to improve precision and calculation > > > stability, please how should I handle the two equations in f and dfdp? > > > > Sorry, please forget what I said about the reshaping of 'x', I was > > confused... I probably shouldn't try to answer help requests while I'm > > working on something different. And I probably don't understand what > > you mean with 'the 2D case'. > > > > The Jacobian should just compute one row of partial derivatives for > > each element of the result of 'f'. > > Thanks but my f returns a 2D matrix (a row for each equation, a column for > each element of the indeps vector x). Should the dfdp return a 3D matrix? From optim_doc, for nonlin_residmin: 'dfdp'      Function computing the Jacobian of the residuals with respect to      the parameters, assuming residuals are reshaped to a column vector. So for calculating the Jacobian, assume that the array returned by 'f' is reshaped to a column vector (as with array(:)). So the Jacobian is 2D. (The differences between nonlin_curvefit and nonlin_residmin do not affect this.) Olaf -- public key id EAFE0591, e.g. on x-hkp://pool.sks-keyservers.net signature.asc (849 bytes) Download Attachment
Open this post in threaded view
|

## Re: Custom jacobian (dfdp) in nonlin_curvefit - calling parameters?

 Dne 18. 11. 19 v 18:39 Olaf Till napsal(a): > > 'dfdp' >       Function computing the Jacobian of the residuals with respect to >       the parameters, assuming residuals are reshaped to a column vector. > > So for calculating the Jacobian, assume that the array returned by 'f' > is reshaped to a column vector (as with array(:)). So the Jacobian is > 2D. > > (The differences between nonlin_curvefit and nonlin_residmin do not > affect this.) > Now my f returns one vector column (first half output of the first function, the second half output of the second function). The default dfdp of minimum differences works perfect. The sum of errors is 1e-17 in about 7 iterations. sum((y - model_values).^2) with 2000 values. My dfdp is IMO correct now, I checked several times. Yet I never get more than 3 iterations, with the sum of errors about 1e-10 which is too inprecise for my case. I configured: settings = optimset ('dfdp', fDfdp, 'TolFun', 1e-17, 'TolX', 1e-17, 'Display', "iter"); No change, always only 2 - 3 iteration. I am listing dfdp output before leaving the function and values in the four columns of dfdp output make sense. Please could the real-finite differences DFDP do someting significantly different than a custom dfdp in terms of settings and tolerances? I did not see anything in __nonlin_residmin__ nor in __lm_svd__ when stepping through the code in debugger. My custom dfdp worked great for values of <0, 1>. But when I am fitting to values < 1e-7 (incl. p values), the custom dfdp is significantly worse than the finite differences. All calculations performed in double. I can stay with the finite differences, but wanted to try if a custom dfdp would further improve the precision since even tiny imprecisions are important in my project. I very much appreciate any recommendation for further troubleshooting. Regards, Pavel.