10:12 p.m.
Dear all, First, a couple of numbers to introduce the question
cos(0.1)~ 0.99500416 5278 025766095561987803870294838576225415084035959352...
### NOTE all results are on 32-bit Ubuntu 15.10 To compare a helper:
function string as_string (scalar x) sprintf s "%.15g",x return s end function
1) x={0.1} eval as_string(fdjac(x,sin(x))) 0.99500416 7780 28
2)
set fdjac_quality 1 eval as_string(fdjac(x,sin(x))) 0.99500416 5451 974
3)
set fdjac_quality 2 eval as_string(fdjac(x,sin(x))) 0.99500416 3123 667
A home-backed thing inspired by fdjac help :
A helper: function scalar feval (string fun, scalar x) return @fun end function
The main:
function scalar sdj2 (string f, scalar x) if abs(x)>10^3 h = 10^3*sqrt(x/1000*$macheps)
else h=10^(3)*sqrt($macheps) endif z = (8*(feval(f,x+h)-feval(f,x-h)) -
(feval(f,x+2*h)-feval(f,x-2*h)))/(12*h) sz = as_string(z) return z end function
We have
eval as_string(sdj2("sin(x)",0.1)) 0.99500416 5277 274 recall: 0.99500416 5278
025766
In theory, fdjac should be better, since it uses the same formula as sdj2() and it should
optimize the step h, but... Oleh
10:48 p.m.
On Tue, 29 Dec 2015, oleg_komashko(a)ukr.net wrote:
Dear all, First, a couple of numbers to introduce the question
cos(0.1)~ 0.99500416 5278 025766095561987803870294838576225415084035959352...
### NOTE all results are on 32-bit Ubuntu 15.10 To compare a helper:
function string as_string (scalar x) sprintf s "%.15g",x return s end function
1) x={0.1} eval as_string(fdjac(x,sin(x))) 0.99500416 7780 28
2)
set fdjac_quality 1 eval as_string(fdjac(x,sin(x))) 0.99500416 5451 974
3)
set fdjac_quality 2 eval as_string(fdjac(x,sin(x))) 0.99500416 3123 667
A home-backed thing inspired by fdjac help :
A helper: function scalar feval (string fun, scalar x) return @fun end function
The main:
function scalar sdj2 (string f, scalar x) if abs(x)>10^3 h =
10^3*sqrt(x/1000*$macheps) else h=10^(3)*sqrt($macheps) endif z =
(8*(feval(f,x+h)-feval(f,x-h)) - (feval(f,x+2*h)-feval(f,x-2*h)))/(12*h) sz = as_string(z)
return z end function
We have
eval as_string(sdj2("sin(x)",0.1)) 0.99500416 5277 274 recall: 0.99500416 5278
025766
In theory, fdjac should be better, since it uses the same formula as sdj2() and it should
optimize the step h, but... Oleh
In general, numerical differentiation should never be expected to be
accurate to more than a few decimal digits, for the combined effect of the
inherent imprecision of machine precision and discrete (as opposed to
_truly_ infinitesimal) differentiation. Run this script and you'll see
what I mean:
<hansl>
set echo off
set messages off
clear
set fdjac_quality 0
loop for (x=-1; x<1; x+=2/13)
matrix a = {x}
printf "x = %5.2f, precision = %d\n", a, \
floor(-log10(abs(1-fdjac(a, sin(a))/cos(x))))
endloop
set fdjac_quality 1
loop for (x=-1; x<1; x+=2/13)
matrix a = {x}
printf "x = %5.2f, precision = %d\n", a, \
floor(-log10(abs(1-fdjac(a, sin(a))/cos(x))))
endloop
set fdjac_quality 2
loop for (x=-1; x<1; x+=2/13)
matrix a = {x}
printf "x = %5.2f, precision = %d\n", a, \
floor(-log10(abs(1-fdjac(a, sin(a))/cos(x))))
endloop
</hansl>
That said, I would have expected Richardson extrapolation to fare somewhat
better in this case. I'll investigate.
-------------------------------------------------------
Riccardo (Jack) Lucchetti
Dipartimento di Scienze Economiche e Sociali (DiSES)
Università Politecnica delle Marche
(formerly known as Università di Ancona)
r.lucchetti(a)univpm.it
http://www2.econ.univpm.it/servizi/hpp/lucchetti
-------------------------------------------------------
11:25 p.m.
Dear Riccardo, My simple script uses exactly Richardson extrapolation but it uses a
simplistic method for defining the step. 3 digits better is by far too large. Clearly,
fdjac uses a more advanced way to find h. It looks, as if it has a buglet inside. I'm
making a wrapper for quadtable to do easy integration with hansl, and my simplistic sfd2*
works systematically better, while it should not to. Oleh *) fd from fdjac
29 грудня 2015, 00:49:01, від "Riccardo (Jack) Lucchetti" <
r.lucchetti(a)univpm.it >:
On Tue, 29 Dec 2015, oleg_komashko(a)ukr.net wrote:
Dear all, First, a couple of numbers to introduce the question
cos(0.1)~ 0.99500416 5278 025766095561987803870294838576225415084035959352...
### NOTE all results are on 32-bit Ubuntu 15.10 To compare a helper:
function string as_string (scalar x) sprintf s "%.15g",x return s end function
1) x={0.1} eval as_string(fdjac(x,sin(x))) 0.99500416 7780 28
2)
set fdjac_quality 1 eval as_string(fdjac(x,sin(x))) 0.99500416 5451 974
3)
set fdjac_quality 2 eval as_string(fdjac(x,sin(x))) 0.99500416 3123 667
A home-backed thing inspired by fdjac help :
A helper: function scalar feval (string fun, scalar x) return @fun end function
The main:
function scalar sdj2 (string f, scalar x) if abs(x)>10^3 h =
10^3*sqrt(x/1000*$macheps) else h=10^(3)*sqrt($macheps) endif z =
(8*(feval(f,x+h)-feval(f,x-h)) - (feval(f,x+2*h)-feval(f,x-2*h)))/(12*h) sz = as_string(z)
return z end function
We have
eval as_string(sdj2("sin(x)",0.1)) 0.99500416 5277 274 recall: 0.99500416 5278
025766
In theory, fdjac should be better, since it uses the same formula as sdj2() and it should
optimize the step h, but... Oleh
In general, numerical differentiation should never be expected to be
accurate to more than a few decimal digits, for the combined effect of the
inherent imprecision of machine precision and discrete (as opposed to
_truly_ infinitesimal) differentiation. Run this script and you'll see
what I mean:
<hansl>
set echo off
set messages off
clear
set fdjac_quality 0
loop for (x=-1; x<1; x+=2/13)
matrix a = {x}
printf "x = %5.2f, precision = %d\n", a, \
floor(-log10(abs(1-fdjac(a, sin(a))/cos(x))))
endloop
set fdjac_quality 1
loop for (x=-1; x<1; x+=2/13)
matrix a = {x}
printf "x = %5.2f, precision = %d\n", a, \
floor(-log10(abs(1-fdjac(a, sin(a))/cos(x))))
endloop
set fdjac_quality 2
loop for (x=-1; x<1; x+=2/13)
matrix a = {x}
printf "x = %5.2f, precision = %d\n", a, \
floor(-log10(abs(1-fdjac(a, sin(a))/cos(x))))
endloop
</hansl>
That said, I would have expected Richardson extrapolation to fare somewhat
better in this case. I'll investigate.
-------------------------------------------------------
Riccardo (Jack) Lucchetti
Dipartimento di Scienze Economiche e Sociali (DiSES)
Università Politecnica delle Marche
(formerly known as Università di Ancona)
r.lucchetti(a)univpm.it
http://www2.econ.univpm.it/servizi/hpp/lucchetti
-------------------------------------------------------
_______________________________________________
Gretl-devel mailing list
Gretl-devel(a)lists.wfu.edu
http://lists.wfu.edu/mailman/listinfo/gretl-devel
1:39 a.m.
On Tue, 29 Dec 2015, oleg_komashko(a)ukr.net wrote:
Dear Riccardo, My simple script uses exactly Richardson
extrapolation but it uses a simplistic method for defining the
step. 3 digits better is by far too large. Clearly, fdjac uses a
more advanced way to find h. It looks, as if it has a buglet
inside.
I don't think there's any bug in Jack's Richardson code. In fact I
believe the only reason you're getting more accurate results in this
case is that you got lucky with your "simplistic method" for setting
h. If I just change the code in our minpack-derived fdjac2.c to use
your choice of
h = 1000 * sqrt(DBL_EPSILON)
in place of
h = fabs(x[j]) * sqrt(DBL_EPSILON)
I get a value for sin'(0.1) which is good to 12 decimal places!
I don't suppose this would be consistent across arbitrary functions,
but it may be worth experimenting with a bigger h, or a different
formula for h, in the "quality = 2" case.
Allin
9:04 a.m.
Dear Allin, For defining h I used the following
if abs(x)>10^3
h = 10^3*abs(x)/1000*sqrt($macheps) else
h=10^(3)*sqrt($macheps) 1000 is more or less recommended by the literature. With
|x|>1000 it will be the same I thought, dfjac uses several calculations to define h.
Interesting. what, say, lapack does use? Oleh
29 грудня 2015, 03:39:57, від "Allin Cottrell" < cottrell(a)wfu.edu >:
On Tue, 29 Dec 2015, oleg_komashko(a)ukr.net wrote:
Dear Riccardo, My simple script uses exactly Richardson
extrapolation but it uses a simplistic method for defining the
step. 3 digits better is by far too large. Clearly, fdjac uses a
more advanced way to find h. It looks, as if it has a buglet
inside.
I don't think there's any bug in Jack's Richardson code. In fact I
believe the only reason you're getting more accurate results in this
case is that you got lucky with your "simplistic method" for setting
h. If I just change the code in our minpack-derived fdjac2.c to use
your choice of
h = 1000 * sqrt(DBL_EPSILON)
in place of
h = fabs(x[j]) * sqrt(DBL_EPSILON)
I get a value for sin'(0.1) which is good to 12 decimal places!
I don't suppose this would be consistent across arbitrary functions,
but it may be worth experimenting with a bigger h, or a different
formula for h, in the "quality = 2" case.
Allin
_______________________________________________
Gretl-devel mailing list
Gretl-devel(a)lists.wfu.edu
http://lists.wfu.edu/mailman/listinfo/gretl-devel
12:14 p.m.
New subject: fdjac: fdjac(x=10^-9,exp(x))=0 instead of 1!
Dear Allin, Thanks for mentioning fdjac2.c: It helped me to find quickly weak spot: Near
zero the step h will be inadequately small
set fdjac_quality 2
x={10^-6}
eval fdjac(x,exp(x))
eval exp(10^-6)
x={10^-7}
eval fdjac(x,exp(x))
eval exp(10^-7)
x={10^-8}
eval fdjac(x,exp(x))
eval exp(10^-8)
x={10^-9}
eval fdjac(x,exp(x))
eval exp(10^-9)
x={10^-12}
eval fdjac(x,exp(x))
eval exp(10^-12)
29 грудня 2015, 03:39:57, від "Allin Cottrell" <cottrell(a)wfu.edu>:
On Tue, 29 Dec 2015, oleg_komashko(a)ukr.net wrote:
> Dear Riccardo, My simple script uses exactly Richardson
> extrapolation but it uses a simplistic method for defining the
> step. 3 digits better is by far too large. Clearly, fdjac uses a
> more advanced way to find h. It looks, as if it has a buglet
> inside.
I don't think there's any bug in Jack's Richardson code. In fact I
believe the only reason you're getting more accurate results in this
case is that you got lucky with your "simplistic method" for setting
h. If I just change the code in our minpack-derived fdjac2.c to use
your choice of
h = 1000 * sqrt(DBL_EPSILON)
in place of
h = fabs(x[j]) * sqrt(DBL_EPSILON)
I get a value for sin'(0.1) which is good to 12 decimal places!
I don't suppose this would be consistent across arbitrary functions,
but it may be worth experimenting with a bigger h, or a different
formula for h, in the "quality = 2" case.
Allin
_______________________________________________
Gretl-devel mailing list
Gretl-devel(a)lists.wfu.edu
http://lists.wfu.edu/mailman/listinfo/gretl-devel
1:34 p.m.
New subject: fdjac: a better example
set fdjac_quality 2 x={10^-9}
function scalar f (scalar x)
return x+10000
end function
eval fdjac(x,f(x))
On 32-bit Ubuntu I have 20345.
Oleh
5:38 p.m.
New subject: fdjac: a better example
On Tue, 29 Dec 2015, oleg_komashko(a)ukr.net wrote:
set fdjac_quality 2
x={10^-9}
function scalar f (scalar x)
return x+10000
end function
eval fdjac(x,f(x))
On 32-bit Ubuntu I have 20345.
Yes, with the status-quo small step size in fdjac2.c, plus
Richardson, the results can become violently unstable for small
values of the function argument. Further investigation:
<hansl>
set echo off
set messages off
function scalar fbig (scalar y)
return y + 10000
end function
scalar f0 f1 f2
loop for (x=1.0; x>1e-10; x/=10) -q
y = {x}
set fdjac_quality 0
f0 = fdjac(y,fbig(y))
set fdjac_quality 1
f1 = fdjac(y,fbig(y))
set fdjac_quality 2
f2 = fdjac(y,fbig(y))
printf "y=%.2e: %f, %f, %f\n", x, f0, f1, f2
endloop
</hansl>
Here, this gives
<gretl-output>
y=1.00e+00: 1.000000, 1.000000, 1.000000
y=1.00e-01: 0.999756, 1.000061, 0.999858
y=1.00e-02: 1.000977, 1.000977, 1.001994
y=1.00e-03: 0.976563, 0.976563, 0.996908
y=1.00e-04: 1.220703, 0.915527, 0.101725
y=1.00e-05: 0.000000, 0.000000, 2.034505
y=1.00e-06: 0.000000, 0.000000, 40.690104
y=1.00e-07: 0.000000, 0.000000, 0.000000
y=1.00e-08: 0.000000, 0.000000, -2034.505208
y=1.00e-09: 0.000000, 0.000000, 20345.052083
y=1.00e-10: 0.000000, 0.000000, 101725.260417
</gretl-output>
It's unsurprising that the "basic" fdjac methods give an answer of
zero when y gets small enough: y + h + 10000 becomes numerically
indistinguishable from y + 10000 when you only have 64 bits per
floating-point number. But the Richardson output goes crazy in a way
that's not very nice at all.
In current git I've modified fdjac2 to use a step size as in your
script (when Richardson is selected). We need to do more testing on
this in case it goes nasty under some conditions, but here's what
I'm now getting from the script above:
<gretl-output>
y=1.00e+00: 1.000000, 1.000000, 1.000000
y=1.00e-01: 0.999756, 1.000061, 1.000000
y=1.00e-02: 1.000977, 1.000977, 1.000000
y=1.00e-03: 0.976563, 0.976563, 1.000000
y=1.00e-04: 1.220703, 0.915527, 1.000000
y=1.00e-05: 0.000000, 0.000000, 1.000000
y=1.00e-06: 0.000000, 0.000000, 1.000000
y=1.00e-07: 0.000000, 0.000000, 1.000000
y=1.00e-08: 0.000000, 0.000000, 1.000000
y=1.00e-09: 0.000000, 0.000000, 1.000000
y=1.00e-10: 0.000000, 0.000000, 1.000000
</gretl-output>
Allin
6:05 p.m.
New subject: fdjac: a better example
On Tue, 29 Dec 2015, Allin Cottrell wrote:
In current git I've modified fdjac2 to use a step size as in your
script
(when Richardson is selected). We need to do more testing on this in case it
goes nasty under some conditions
It would seem that the results are somewhat sensitive to the amount of
nonlinearity of the problem: I ran an old script I had which uses matrix
inversion (notoriously differentiable almost everywhere but highly
nonlinear --- especially in the nieighbourhood of singularity).
<hansl>
set echo off
set messages off
function matrix vec22inv(matrix x)
n = round(sqrt(rows(x)))
A = mshape(x,n,n)
return vec(inv(A))'
end function
function matrix grad(matrix x)
n = round(sqrt(rows(x)))
iA = inv(mshape(x,n,n))
return -iA' ** iA
end function
n = 12
k = 1.0e4
a = vec(ones(n,n) - I(n)/k)
ia = vec22inv(a)
ag = grad(a)
printf "\n"
loop i=0..2 --quiet
set fdjac_quality $i
set stopwatch
ng$i = fdjac(a, "vec22inv(a)")
tt = $stopwatch
scalar aae = meanc(abs(vec(ag - ng$i)))
scalar mae = maxc(abs(vec(ag - ng$i)))
scalar are = meanc(abs(1 - vec(ag ./ ng$i)))
scalar mre = maxc(abs(1 - vec(ag ./ ng$i)))
printf "quality = $i\n"
printf "\tAverage and max absolute error = %g, %g\n", aae, mae
printf "\tAverage and max relative error = %g, %g\n", are, mre
printf "\t(cpu time = %g)\n", tt
endloop
</hansl>
This is what I get before and after the patch:
<before>
quality = 0
Average and max absolute error = 53.1382, 11477.9
Average and max relative error = 2.27644e-05, 0.00013658
(cpu time = 0.02)
quality = 1
Average and max absolute error = 0.0163915, 6.19778
Average and max relative error = 7.58353e-09, 1.68055e-07
(cpu time = 0.04)
quality = 2
Average and max absolute error = 0.00128204, 0.149397
Average and max relative error = 1.11463e-09, 2.15128e-07
(cpu time = 0.09)
</before>
<after>
quality = 0
Average and max absolute error = 53.1382, 11477.9
Average and max relative error = 2.27644e-05, 0.00013658
(cpu time = 0.02)
quality = 1
Average and max absolute error = 0.0163915, 6.19778
Average and max relative error = 7.58353e-09, 1.68055e-07
(cpu time = 0.08)
quality = 2
Average and max absolute error = 298.462, 128846
Average and max relative error = 0.000128065, 0.00153573
(cpu time = 0.07)
</after>
-------------------------------------------------------
Riccardo (Jack) Lucchetti
Dipartimento di Scienze Economiche e Sociali (DiSES)
Università Politecnica delle Marche
(formerly known as Università di Ancona)
r.lucchetti(a)univpm.it
http://www2.econ.univpm.it/servizi/hpp/lucchetti
-------------------------------------------------------
1:08 a.m.
New subject: fdjac: a better example
On Tue, 29 Dec 2015, Riccardo (Jack) Lucchetti wrote:
On Tue, 29 Dec 2015, Allin Cottrell wrote:
> In current git I've modified fdjac2 to use a step size as in your script
> (when Richardson is selected). We need to do more testing on this in case
> it goes nasty under some conditions
It would seem that the results are somewhat sensitive to the amount of
nonlinearity of the problem: I ran an old script I had which uses matrix
inversion (notoriously differentiable almost everywhere but highly nonlinear
--- especially in the nieighbourhood of singularity).
<hansl>
set echo off
set messages off
function matrix vec22inv(matrix x)
n = round(sqrt(rows(x)))
A = mshape(x,n,n)
return vec(inv(A))'
end function
function matrix grad(matrix x)
n = round(sqrt(rows(x)))
iA = inv(mshape(x,n,n))
return -iA' ** iA
end function
n = 12
k = 1.0e4
a = vec(ones(n,n) - I(n)/k)
ia = vec22inv(a)
ag = grad(a)
printf "\n"
loop i=0..2 --quiet
set fdjac_quality $i
set stopwatch
ng$i = fdjac(a, "vec22inv(a)")
tt = $stopwatch
scalar aae = meanc(abs(vec(ag - ng$i)))
scalar mae = maxc(abs(vec(ag - ng$i)))
scalar are = meanc(abs(1 - vec(ag ./ ng$i)))
scalar mre = maxc(abs(1 - vec(ag ./ ng$i)))
printf "quality = $i\n"
printf "\tAverage and max absolute error = %g, %g\n", aae, mae
printf "\tAverage and max relative error = %g, %g\n", are, mre
printf "\t(cpu time = %g)\n", tt
endloop
</hansl>
This is what I get before and after the patch:
<before>
quality = 0
Average and max absolute error = 53.1382, 11477.9
Average and max relative error = 2.27644e-05, 0.00013658
(cpu time = 0.02)
quality = 1
Average and max absolute error = 0.0163915, 6.19778
Average and max relative error = 7.58353e-09, 1.68055e-07
(cpu time = 0.04)
quality = 2
Average and max absolute error = 0.00128204, 0.149397
Average and max relative error = 1.11463e-09, 2.15128e-07
(cpu time = 0.09)
</before>
<after>
quality = 0
Average and max absolute error = 53.1382, 11477.9
Average and max relative error = 2.27644e-05, 0.00013658
(cpu time = 0.02)
quality = 1
Average and max absolute error = 0.0163915, 6.19778
Average and max relative error = 7.58353e-09, 1.68055e-07
(cpu time = 0.08)
quality = 2
Average and max absolute error = 298.462, 128846
Average and max relative error = 0.000128065, 0.00153573
(cpu time = 0.07)
</after>
Interesting. The "after" case is pretty compelling evidence that
Oleh's setting of the step size -- which worked so nicely in some
cases -- can produce very bad results in others.
However, here (building with gcc 5.3 on Fedora 23), the results are
not so impressive for the "quality 2" case under the status quo ante
(git commit ee7808 of fdjac2.c). I'm seeing:
<before>
quality = 0
Average and max absolute error = 53.1431, 11477.9
Average and max relative error = 2.27665e-05, 0.000136579
(cpu time = 0.0192254)
quality = 1
Average and max absolute error = 0.0448263, 6.1979
Average and max relative error = 3.92204e-08, 8.20507e-07
(cpu time = 0.0307975)
quality = 2
Average and max absolute error = 0.0814845, 1.41264
Average and max relative error = 9.3409e-08, 2.03417e-06
(cpu time = 0.073777)
</before>
Do we have a case where "quality 2" unambiguously produces better
results than "quality 1"?
Allin
8:36 a.m.
New subject: fdjac: a better example
We have 2 (3) ways of proceed : make experiments, find good expert on numerical methods
(do both) The script attached shows, for small x the best h is as big as
10^5*sqrt($macheps) Still I do not know what to to with x: abs(x)>10^5 and what to do
with say, exp(100) try x = {100} eval fdjac (x,exp(x)) - exp(100) That's ERROR! Oleh
30 грудня 2015, 03:08:30, від "Allin Cottrell" < cottrell(a)wfu.edu >:
On Tue, 29 Dec 2015, Riccardo (Jack) Lucchetti wrote:
On Tue, 29 Dec 2015, Allin Cottrell wrote:
> In current git I've modified fdjac2 to use a step size as in your script
> (when Richardson is selected). We need to do more testing on this in case
> it goes nasty under some conditions
It would seem that the results are somewhat sensitive to the amount of
nonlinearity of the problem: I ran an old script I had which uses matrix
inversion (notoriously differentiable almost everywhere but highly nonlinear
--- especially in the nieighbourhood of singularity).
<hansl>
set echo off
set messages off
function matrix vec22inv(matrix x)
n = round(sqrt(rows(x)))
A = mshape(x,n,n)
return vec(inv(A))'
end function
function matrix grad(matrix x)
n = round(sqrt(rows(x)))
iA = inv(mshape(x,n,n))
return -iA' ** iA
end function
n = 12
k = 1.0e4
a = vec(ones(n,n) - I(n)/k)
ia = vec22inv(a)
ag = grad(a)
printf "\n"
loop i=0..2 --quiet
set fdjac_quality $i
set stopwatch
ng$i = fdjac(a, "vec22inv(a)")
tt = $stopwatch
scalar aae = meanc(abs(vec(ag - ng$i)))
scalar mae = maxc(abs(vec(ag - ng$i)))
scalar are = meanc(abs(1 - vec(ag ./ ng$i)))
scalar mre = maxc(abs(1 - vec(ag ./ ng$i)))
printf "quality = $i\n"
printf "\tAverage and max absolute error = %g, %g\n", aae, mae
printf "\tAverage and max relative error = %g, %g\n", are, mre
printf "\t(cpu time = %g)\n", tt
endloop
</hansl>
This is what I get before and after the patch:
<before>
quality = 0
Average and max absolute error = 53.1382, 11477.9
Average and max relative error = 2.27644e-05, 0.00013658
(cpu time = 0.02)
quality = 1
Average and max absolute error = 0.0163915, 6.19778
Average and max relative error = 7.58353e-09, 1.68055e-07
(cpu time = 0.04)
quality = 2
Average and max absolute error = 0.00128204, 0.149397
Average and max relative error = 1.11463e-09, 2.15128e-07
(cpu time = 0.09)
</before>
<after>
quality = 0
Average and max absolute error = 53.1382, 11477.9
Average and max relative error = 2.27644e-05, 0.00013658
(cpu time = 0.02)
quality = 1
Average and max absolute error = 0.0163915, 6.19778
Average and max relative error = 7.58353e-09, 1.68055e-07
(cpu time = 0.08)
quality = 2
Average and max absolute error = 298.462, 128846
Average and max relative error = 0.000128065, 0.00153573
(cpu time = 0.07)
</after>
Interesting. The "after" case is pretty compelling evidence that
Oleh's setting of the step size -- which worked so nicely in some
cases -- can produce very bad results in others.
However, here (building with gcc 5.3 on Fedora 23), the results are
not so impressive for the "quality 2" case under the status quo ante
(git commit ee7808 of fdjac2.c). I'm seeing:
<before>
quality = 0
Average and max absolute error = 53.1431, 11477.9
Average and max relative error = 2.27665e-05, 0.000136579
(cpu time = 0.0192254)
quality = 1
Average and max absolute error = 0.0448263, 6.1979
Average and max relative error = 3.92204e-08, 8.20507e-07
(cpu time = 0.0307975)
quality = 2
Average and max absolute error = 0.0814845, 1.41264
Average and max relative error = 9.3409e-08, 2.03417e-06
(cpu time = 0.073777)
</before>
Do we have a case where "quality 2" unambiguously produces better
results than "quality 1"?
Allin
_______________________________________________
Gretl-devel mailing list
Gretl-devel(a)lists.wfu.edu
http://lists.wfu.edu/mailman/listinfo/gretl-devel
10:48 a.m.
New subject: fdjac: a better example
On Wed, 30 Dec 2015, oleg_komashko(a)ukr.net wrote:
We have 2 (3) ways of proceed : make experiments, find good expert on
numerical methods (do both) The script attached shows, for small x the
best h is as big as 10^5*sqrt($macheps) Still I do not know what to to
with x: abs(x)>10^5 and what to do with say, exp(100) try x = {100} eval
fdjac (x,exp(x)) - exp(100) That's ERROR! Oleh
Of course ANY numerical approximation to a derivative is a ratio;
therefore, what is crucial to its precision is the relative order of
magnitude of numerator and denominator. It's easy to contruct pathological
cases and there's no silver bullet that works in all cases.
That said, the following settings seem to work reasonably well in a
decently-sized array of cases (diff against git HEAD --- not applied yet).
<diff>
diff --git a/minpack/fdjac2.c b/minpack/fdjac2.c
index 9e9676a..54e64ec 100644
--- a/minpack/fdjac2.c
+++ b/minpack/fdjac2.c
@@ -144,10 +144,10 @@ int fdjac2_(S_fp fcn, int m, int n, int quality,
for (j = 0; j < n; j++) {
temp = x[j];
if (quality == 2) {
- if (fabs(temp) > 1000) {
- h = 1000 * sqrt(fabs(temp/1000)) * eps;
+ if (fabs(temp) > 128) {
+ h = 128* sqrt(fabs(temp/128)) * eps / sqrt(6);
} else {
- h = 1000 * eps;
+ h = 128 * eps / sqrt(6);
}
} else {
h = eps * fabs(temp);
</diff>
I'm attaching a battery of tests.
-------------------------------------------------------
Riccardo (Jack) Lucchetti
Dipartimento di Scienze Economiche e Sociali (DiSES)
Università Politecnica delle Marche
(formerly known as Università di Ancona)
r.lucchetti(a)univpm.it
http://www2.econ.univpm.it/servizi/hpp/lucchetti
-------------------------------------------------------
11:48 a.m.
New subject: fdjac: a better example
Of course, my exp(100) example is incorrect: relative precision is better than 10^-12 even
on 32-bit system. Also, ways to determine h should be different for qualities 0,1, and 2
Quality 0 requires much lesser h: for x from (0,1) sqrt($macheps) even 10^-1*
sqrt($macheps) works OK And of course if x==0 should be changed to something like
abs(x)< 10^-1* sqrt($macheps) for quality 0 Currently, I can't tell about quality 1
Oleh I'll write my impression on the tests later
30 грудня 2015, 12:48:19, від "Riccardo (Jack) Lucchetti" <
r.lucchetti(a)univpm.it >:
On Wed, 30 Dec 2015, oleg_komashko(a)ukr.net wrote:
We have 2 (3) ways of proceed : make experiments, find good expert on
numerical methods (do both) The script attached shows, for small x the
best h is as big as 10^5*sqrt($macheps) Still I do not know what to to
with x: abs(x)>10^5 and what to do with say, exp(100) try x = {100} eval
fdjac (x,exp(x)) - exp(100) That's ERROR! Oleh
Of course ANY numerical approximation to a derivative is a ratio;
therefore, what is crucial to its precision is the relative order of
magnitude of numerator and denominator. It's easy to contruct pathological
cases and there's no silver bullet that works in all cases.
That said, the following settings seem to work reasonably well in a
decently-sized array of cases (diff against git HEAD --- not applied yet).
<diff>
diff --git a/minpack/fdjac2.c b/minpack/fdjac2.c
index 9e9676a..54e64ec 100644
--- a/minpack/fdjac2.c
+++ b/minpack/fdjac2.c
@@ -144,10 +144,10 @@ int fdjac2_(S_fp fcn, int m, int n, int quality,
for (j = 0; j < n; j++) {
temp = x[j];
if (quality == 2) {
- if (fabs(temp) > 1000) {
- h = 1000 * sqrt(fabs(temp/1000)) * eps;
+ if (fabs(temp) > 128) {
+ h = 128* sqrt(fabs(temp/128)) * eps / sqrt(6);
} else {
- h = 1000 * eps;
+ h = 128 * eps / sqrt(6);
}
} else {
h = eps * fabs(temp);
</diff>
I'm attaching a battery of tests.
-------------------------------------------------------
Riccardo (Jack) Lucchetti
Dipartimento di Scienze Economiche e Sociali (DiSES)
Università Politecnica delle Marche
(formerly known as Università di Ancona)
r.lucchetti(a)univpm.it
http://www2.econ.univpm.it/servizi/hpp/lucchetti
-------------------------------------------------------
2:29 p.m.
New subject: fdjac: a better example
On tests: look reasonably well
For now I see 2 points: 1) Since 10^3*sqrt($macheps) is ~10^-5 on 32 bit systems there is
a problem with, e.g. sqrt(x) near zero. May be, simply make warning in the help text? 2)
strictly speaking sqrt(abs(x)) is incorrect and should be replaced by a linear function of
abs(x). Still, I don't know what exactly it should be. For quality 0 the following
works reasonably well: (for sin(x) for x in (0, 5*10^7) ) if abs(x)<1
h=sqrt($macheps)
else h= (1+(2^(-12))*(abs(x)-1))*sqrt($macheps) endif
Also, I think it will be nice to also have the left-sided version of quality 0: for 1/x
with small negative x, for example Oleh
30 грудня 2015, 12:48:19, від "Riccardo (Jack) Lucchetti" <
r.lucchetti(a)univpm.it >:
On Wed, 30 Dec 2015, oleg_komashko(a)ukr.net wrote:
We have 2 (3) ways of proceed : make experiments, find good expert on
numerical methods (do both) The script attached shows, for small x the
best h is as big as 10^5*sqrt($macheps) Still I do not know what to to
with x: abs(x)>10^5 and what to do with say, exp(100) try x = {100} eval
fdjac (x,exp(x)) - exp(100) That's ERROR! Oleh
Of course ANY numerical approximation to a derivative is a ratio;
therefore, what is crucial to its precision is the relative order of
magnitude of numerator and denominator. It's easy to contruct pathological
cases and there's no silver bullet that works in all cases.
That said, the following settings seem to work reasonably well in a
decently-sized array of cases (diff against git HEAD --- not applied yet).
<diff>
diff --git a/minpack/fdjac2.c b/minpack/fdjac2.c
index 9e9676a..54e64ec 100644
--- a/minpack/fdjac2.c
+++ b/minpack/fdjac2.c
@@ -144,10 +144,10 @@ int fdjac2_(S_fp fcn, int m, int n, int quality,
for (j = 0; j < n; j++) {
temp = x[j];
if (quality == 2) {
- if (fabs(temp) > 1000) {
- h = 1000 * sqrt(fabs(temp/1000)) * eps;
+ if (fabs(temp) > 128) {
+ h = 128* sqrt(fabs(temp/128)) * eps / sqrt(6);
} else {
- h = 1000 * eps;
+ h = 128 * eps / sqrt(6);
}
} else {
h = eps * fabs(temp);
</diff>
I'm attaching a battery of tests.
-------------------------------------------------------
Riccardo (Jack) Lucchetti
Dipartimento di Scienze Economiche e Sociali (DiSES)
Università Politecnica delle Marche
(formerly known as Università di Ancona)
r.lucchetti(a)univpm.it
http://www2.econ.univpm.it/servizi/hpp/lucchetti
-------------------------------------------------------
3:01 p.m.
New subject: fdjac: a better example
On Wed, 30 Dec 2015, oleg_komashko(a)ukr.net wrote:
Also, I think it will be nice to also have the left-sided version of
quality 0: for 1/x with small negative x, for example
really? can you make an example in which it would be useful? (apart from
teaching purposes?)
-------------------------------------------------------
Riccardo (Jack) Lucchetti
Dipartimento di Scienze Economiche e Sociali (DiSES)
Università Politecnica delle Marche
(formerly known as Università di Ancona)
r.lucchetti(a)univpm.it
http://www2.econ.univpm.it/servizi/hpp/lucchetti
-------------------------------------------------------
3:10 p.m.
New subject: fdjac: a better example
OK, I do not insist But it's more interesting, what do you think on changing h for 0
and 1 qualities?
30 грудня 2015, 17:01:34, від "Riccardo (Jack) Lucchetti" <
r.lucchetti(a)univpm.it >:
On Wed, 30 Dec 2015, oleg_komashko(a)ukr.net wrote:
Also, I think it will be nice to also have the left-sided version of
quality 0: for 1/x with small negative x, for example
really? can you make an example in which it would be useful? (apart from
teaching purposes?)
-------------------------------------------------------
Riccardo (Jack) Lucchetti
Dipartimento di Scienze Economiche e Sociali (DiSES)
Università Politecnica delle Marche
(formerly known as Università di Ancona)
r.lucchetti(a)univpm.it
http://www2.econ.univpm.it/servizi/hpp/lucchetti
-------------------------------------------------------
6:50 p.m.
New subject: fdjac: some quality 0 tests
1) Thank you for BIGGER_H 2) OK Now I know where tiny abs(x)*sqrt($macheps) for x's
near zero are better: x^alpha, alpha<1 Now I can use such way around: use quality 1 for
powers<1 near zero and quality 2 for the rest Oleh P.S. I was surprised to see such
established thing as Minpack does the way it does. I think, the legacy code must contain
something better. Of course, the problem is closed for the time being.
30 грудня 2015, 20:38:25, від "Allin Cottrell" < cottrell(a)wfu.edu >:
On Wed, 30 Dec 2015, oleg_komashko(a)ukr.net wrote:
Attached
Our "quality 0" is the standard Minpack code, I don't think we want
to mess with it. (But we may decide it shouldn't be the default.)
Allin
2:44 p.m.
New subject: fdjac: a better example
On Wed, 30 Dec 2015, Riccardo (Jack) Lucchetti wrote:
On Wed, 30 Dec 2015, oleg_komashko(a)ukr.net wrote:
> We have 2 (3) ways of proceed : make experiments, find good expert on
> numerical methods (do both) The script attached shows, for small x the best
> h is as big as 10^5*sqrt($macheps) Still I do not know what to to with x:
> abs(x)>10^5 and what to do with say, exp(100) try x = {100} eval fdjac
> (x,exp(x)) - exp(100) That's ERROR! Oleh
Of course ANY numerical approximation to a derivative is a ratio; therefore,
what is crucial to its precision is the relative order of magnitude of
numerator and denominator. It's easy to contruct pathological cases and
there's no silver bullet that works in all cases.
That said, the following settings seem to work reasonably well in a
decently-sized array of cases (diff against git HEAD --- not applied yet).
<diff>
diff --git a/minpack/fdjac2.c b/minpack/fdjac2.c
index 9e9676a..54e64ec 100644
--- a/minpack/fdjac2.c
+++ b/minpack/fdjac2.c
@@ -144,10 +144,10 @@ int fdjac2_(S_fp fcn, int m, int n, int quality,
for (j = 0; j < n; j++) {
temp = x[j];
if (quality == 2) {
- if (fabs(temp) > 1000) {
- h = 1000 * sqrt(fabs(temp/1000)) * eps;
+ if (fabs(temp) > 128) {
+ h = 128* sqrt(fabs(temp/128)) * eps / sqrt(6);
} else {
- h = 1000 * eps;
+ h = 128 * eps / sqrt(6);
}
} else {
h = eps * fabs(temp);
</diff>
I'm attaching a battery of tests.
Thanks, Jack. This does look like a reasonable compromise. For ease
of experimentation I've committed a version of your patch, governed
by the BIGGER_H define (set on for the moment).
Allin
6:27 p.m.
New subject: fdjac: a better example
On Wed, 30 Dec 2015, oleg_komashko(a)ukr.net wrote:
> governed by the BIGGER_H define (set on for the moment) Excuse
> me, this is about c code, or I have to something more than set
> fdjac_quality 2?
The BIGGER_H thing is in fdjac2.c. Right now it's set to give the
variant that Jack posted recently when fdjac_quality 2 is selected.
Allin
3280
days inactive
3282
days old
21 comments
3 participants
participants (3)
-
Allin Cottrell -
oleg_komashko＠ukr.net -
Riccardo (Jack) Lucchetti