

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);


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 mcode (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 xhkp://pool.skskeyservers.net


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 mcode (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.


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 mcode (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
elementbyelement 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.


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 xhkp://pool.skskeyservers.net


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.


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 xhkp://pool.skskeyservers.net


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.


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 xhkp://pool.skskeyservers.net


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 1e17 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 1e10 which is too
inprecise for my case.
I configured:
settings = optimset ('dfdp', fDfdp, 'TolFun', 1e17, 'TolX', 1e17,
'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 realfinite 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 < 1e7 (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.


On Tue, Nov 19, 2019 at 12:00:25AM +0100, Pavel Hofman wrote:
>
> 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
> 1e17 in about 7 iterations.
>
> sum((y  model_values).^2) with 2000 values.
This tells us nothing without knowing how large your elements of 'y'
are.
> My dfdp is IMO correct now, I checked several times. Yet I never get more
> than 3 iterations, with the sum of errors about 1e10 which is too inprecise
> for my case.
My guess is that your custom dfdp still has some error (I didn't
check).
You could check your custom dfdp against finite differences, maybe
playing with 'FinDiffRelStep', 'FinDiffType', and 'TypicalX'. For more
accuracy you could use the function 'jacobs' for complex step
derivatives, if your 'f' allows it.
A study of 'optim_doc' should help you in figuring out the above.
> I configured:
>
> settings = optimset ('dfdp', fDfdp, 'TolFun', 1e17, 'TolX', 1e17,
> 'Display', "iter");
'TolX' is not honoured by the default algorithm (which is mentioned by
'optim_doc'). 'TolFun' is the minimum _fractional_ improvement of the
sum of squared residuals, so setting it to less than machine epsilon,
as you did, means the algorithm will start a new outer loop after any
improvement (which is not necessarily good, since the inner loop won't
further try to refine the direction of the step).
> 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 realfinite 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.
No, the only thing different should be which Jacobian function is
called.
Olaf

public key id EAFE0591, e.g. on xhkp://pool.skskeyservers.net


Hi Olaf,
Dne 19. 11. 19 v 12:04 Olaf Till napsal(a):
>
> My guess is that your custom dfdp still has some error (I didn't
> check).
Thanks a lot for your patience. Of course I had wrong order of p
elements in my dfdp. The custom precisely calculated jacobian yields
better results than the finite differences method now, consistent with
all the other cases in my project.
Once more thank you.
Best regards,
Pavel.

