[ADMB Users] "Hessian does not appear to be positive definite"
Sylvain Bonhommeau
sylvain.bonhommeau at gmail.com
Thu Apr 9 01:20:34 PDT 2009
Hi users,
I played with ADMB to test it before my real application. I first simulated
data resulting from the model I would like to test. It is a weird non-linear
model with 24 parameters (so 25 with sigma). JFY, the general expression for
this model (in R):
Y[i]<-A[i]+(w0+w1*X+w2*X^2)*exp(-(s0+s1*X+s2*X^2)*(B[i]-443))+a1[i]*V^a2[i]+b*Z*(B[i]/443)^n
(don't worry it's an ocean optic model...)
"i" is the wavelength of light and there are 8 wavelengths
A[i], B[i] are constant depending on i
X, V and Z are my independent variables
and w0, w1, w2, s0, s1, s2, a1[i], a2[i], b, n are the parameters
I fixed the parameters and simulate Y to test the ADMB script. The retrieved
parameters have exactly (and amazingly!) the same values as I fixed (which
is not the case when using optim function in R or fminsearch in matlab and
it's much longer).
The problem is that I have the "error" message: "Hessian does not appear to
be positive definite" and I cannot get a "sdreport" for the parameters or
the variable Y.
I copy paste the last part of the ADMB run as well as the tpl file
hereafter.
Thanks for your help,
Sylvain
25 variables; iteration 740; function evaluation 907
Function value -2.7878808e+05; maximum gradient component mag 5.7352e+18
Var Value Gradient |Var Value Gradient |Var Value
Gradient
1-36.23879 4.58537e+02 | 2 0.00190 -3.50573e+18 | 3 1.30460
2.70468e+16
4 1.37629 7.65838e+15 | 5 0.01472 5.73518e+18 | 6 0.00288
1.15859e+18
7 -0.03205 1.47477e+17 | 8 -0.52865 -4.26941e+14 | 9 -0.52865
1.04332e+15
10 -1.47135 -6.57274e+14 | 11 0.12112 -1.61428e+16 | 12 0.16213
-2.80633e+16
13 0.08573 -1.75030e+16 | 14 -0.00369 -1.96681e+16 | 15 -0.21428
7.12243e+15
16-11.81483 1.82708e+14 | 17 4.18517 2.72772e+14 | 18 1.81483
3.84233e+13
19 0.21166 3.25505e+13 | 20 0.09244 1.03861e+16 | 21 0.05250
1.08502e+15
22 0.26054 2.77481e+15 | 23 0.78956 -9.50897e+13 | 24 210.1177
2.82176e+09
25 227.8695 -9.76063e+06 |
ic > imax in fminim is answer attained ?
Var Value Gradient |Var Value Gradient |Var Value
Gradient
1-36.23882 4.58537e+02 | 2 0.00190 -3.50573e+18 | 3 1.30460
2.70468e+16
4 1.37629 7.65838e+15 | 5 0.01472 5.73518e+18 | 6 0.00288
1.15859e+18
7 -0.03205 1.47477e+17 | 8 -0.52865 -4.26941e+14 | 9 -0.52865
1.04332e+15
10 -1.47135 -6.57274e+14 | 11 0.12112 -1.61428e+16 | 12 0.16213
-2.80633e+16
13 0.08573 -1.75030e+16 | 14 -0.00369 -1.96681e+16 | 15 -0.21428
7.12243e+15
16-11.81483 1.82708e+14 | 17 4.18517 2.72772e+14 | 18 1.81483
3.84233e+13
19 0.21166 3.25505e+13 | 20 0.09244 1.03861e+16 | 21 0.05250
1.08502e+15
22 0.26054 2.77481e+15 | 23 0.78956 -9.50897e+13 | 24 210.1177
2.82176e+09
25 227.8695 -9.76063e+06 |
Function minimizer not making progress ... is minimum attained?
Minimprove criterion = 0.0000e+00
- final statistics:
25 variables; iteration 741; function evaluation 937
Function value -2.7879e+05; maximum gradient component mag 5.7352e+18
Exit code = 1; converg criter 1.0000e-12
Var Value Gradient |Var Value Gradient |Var Value
Gradient
1-36.23882 4.58537e+02 | 2 0.00190 -3.50573e+18 | 3 1.30460
2.70468e+16
4 1.37629 7.65838e+15 | 5 0.01472 5.73518e+18 | 6 0.00288
1.15859e+18
7 -0.03205 1.47477e+17 | 8 -0.52865 -4.26941e+14 | 9 -0.52865
1.04332e+15
10 -1.47135 -6.57274e+14 | 11 0.12112 -1.61428e+16 | 12 0.16213
-2.80633e+16
13 0.08573 -1.75030e+16 | 14 -0.00369 -1.96681e+16 | 15 -0.21428
7.12243e+15
16-11.81483 1.82708e+14 | 17 4.18517 2.72772e+14 | 18 1.81483
3.84233e+13
19 0.21166 3.25505e+13 | 20 0.09244 1.03861e+16 | 21 0.05250
1.08502e+15
22 0.26054 2.77481e+15 | 23 0.78956 -9.50897e+13 | 24 210.1177
2.82176e+09
25 227.8695 -9.76063e+06 |
Estimating row 1 out of 25 for hessian
Estimating row 2 out of 25 for hessian
Estimating row 3 out of 25 for hessian
Estimating row 4 out of 25 for hessian
Estimating row 5 out of 25 for hessian
Estimating row 6 out of 25 for hessian
Estimating row 7 out of 25 for hessian
Estimating row 8 out of 25 for hessian
Estimating row 9 out of 25 for hessian
Estimating row 10 out of 25 for hessian
Estimating row 11 out of 25 for hessian
Estimating row 12 out of 25 for hessian
Estimating row 13 out of 25 for hessian
Estimating row 14 out of 25 for hessian
Estimating row 15 out of 25 for hessian
Estimating row 16 out of 25 for hessian
Estimating row 17 out of 25 for hessian
Estimating row 18 out of 25 for hessian
Estimating row 19 out of 25 for hessian
Estimating row 20 out of 25 for hessian
Estimating row 21 out of 25 for hessian
Estimating row 22 out of 25 for hessian
Estimating row 23 out of 25 for hessian
Estimating row 24 out of 25 for hessian
Estimating row 25 out of 25 for hessian
Warning -- Hessian does not appear to be positive definite
Hessian does not appear to be positive definite
///////////////////////////////////////////////////////////////////////////////////////////////////////////
TPLfile
///////////////////////////////////////////////////////////////////////////////////////////////////////////
DATA_SECTION
init_int Nobs_tot //Total number of observations
init_vector Kw(1,8) //Value of Water absorption for each
wavelength
init_matrix Data(1,Nobs_tot,1,6) // Data arranged in column such as
:Kd, Chl, Acdm, Bbp, lambda, convert_lambda
vector Kd(1,Nobs_tot) // Kd data as a vector
vector Chl(1,Nobs_tot) // Chlorophyll data as a vector
vector Acdm(1,Nobs_tot) // Acdm data as a vector
vector Bbp(1,Nobs_tot) // Bbp data as a vector
vector lambda(1,Nobs_tot) // Lambda as vector
vector convert_lambda(1,Nobs_tot) // Lambda as an indice 1=320, 2=340,
3=380, 4=412, 5=443, 6=490, 7=510, 8=555
PARAMETER_SECTION
init_number logSigma;
number sigma;
init_vector w(1,3) //3 common parameters for the
quadratic function of acdm (equivalent to CC)
init_vector s(1,3) //3 common parameters for the
quadratic function of acdm in the exponent term (equivalent to S)
init_bounded_vector loga1(1,8,-7,-1) //8 wavelength-specific
parameters that HAVE to be positive
init_bounded_vector loga2(1,8,-1,-0.0001) //8
wavelength-specific parameters (as an exponent) that Have to be positive
init_bounded_number logbeta0(-20,-10) //1 common parameter
for the Bbp that HAS to be positive
init_bounded_number logeta(-0.7,1.1) //1 common parameter
for the power of wavelength ratio that HAS to be positive
vector pred_Kd(1,Nobs_tot) //one vector of the predicted value
of Kd
vector a1(1,8)
sdreport_vector a3(1,8)
vector a2(1,8)
number beta0
number eta
objective_function_value f
PRELIMINARY_CALCS_SECTION
Kd=column(Data,1);
Chl=column(Data,2);
Acdm=column(Data,3);
Bbp=column(Data,4);
lambda=column(Data,5);
convert_lambda=column(Data,6);
PROCEDURE_SECTION
f=0;
sigma=exp(logSigma);
a1=exp(loga1);
a2=exp(loga2);
a3=a1;
beta0=exp(logbeta0);
eta=-exp(logeta);
for (int i=1;i<=Nobs_tot;i++)
{
pred_Kd(i)=Kw(convert_lambda(i))+(w(1)+w(2)*Acdm(i)+w(3)*Acdm(i)*Acdm(i))*exp(-(s(1)+s(2)*Acdm(i)+s(3)*Acdm(i)*Acdm(i))*(lambda(i)-443.0))+a1(convert_lambda(i))*pow(Chl(i),a2(convert_lambda(i)))+beta0*Bbp(i)*pow((lambda(i)/443.0),eta);
}
// cout<<sum(square(pred_Kd-Kd))<<endl;
f=0.5*(Nobs_tot*log(2*M_PI*square(sigma))+sum(square(pred_Kd-Kd)/square(sigma)));
RUNTIME_SECTION
maximum_function_evaluations 200000000
convergence_criteria 1.e-12
TOP_OF_MAIN_SECTION
arrmblsize = 100000000;
gradient_structure::set_MAX_NVAR_OFFSET(3000);
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.admb-project.org/pipermail/users/attachments/20090409/7aa8d7b0/attachment.html>
More information about the Users
mailing list