[ADMB Users] "Hessian does not appear to be positive definite"

473 views
Skip to first unread message

Sylvain Bonhommeau

unread,
Apr 9, 2009, 4:20:34 AM4/9/09
to us...@admb-project.org
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);

dave fournier

unread,
Apr 9, 2009, 12:08:31 PM4/9/09
to us...@admb-project.org
Hi,

I think there is a bit lost in translation here.

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

Do you mean the you used known parameters to simulate some data and then
fit the data and got the parameters back. Is that the same data set
for which the results you posted here are from? anyway if you send me
the tpl file and every thing needed to run the model I'll take a look at it.

Dave


--
David A. Fournier
P.O. Box 2040,
Sidney, B.C. V8l 3S3
Canada
Phone/FAX 250-655-3364
http://otter-rsch.com
_______________________________________________
Users mailing list
Us...@admb-project.org
http://lists.admb-project.org/mailman/listinfo/users

dave fournier

unread,
Apr 9, 2009, 1:21:32 PM4/9/09
to us...@admb-project.org

Your problme is that the data are too good so the model fits perfectly
with an estimated sigma of almost 0. I added a tiny bit of noise to the
data and it converged fine.

Kd=column(Data,1);
random_number_generator rng(301);
dvector tmp(1,Nobs_tot);
Kd=elem_prod(Kd,exp(.01*tmp)); // add some lognormal noise to Kd

However with real data you may need to civilize the code a bit.

dave fournier

unread,
Apr 9, 2009, 3:24:28 PM4/9/09
to us...@admb-project.org
I think we should keep this on the list as these are issures whihc may
interest other users. If the data are perfect then the residuals are
amost zero and the estimate for the standard deviation is almost zero so
you get a

0/0

situation. very unlikely that this will occur with real data.

If you are moving on to random effects you should first civilize the
model a bit. In somethng like

exp( a(1) + a(2)*v + a(3)*v^2 )

you should at least center v as in

w = v -mean(v)

mfexp( a(1) + a(2)*w + a(3)*w^2 )

The use of mfexp is to avoid overflow.

also you might fit a(1) first as in

PARAMETER_SECTION

init_number a1 // fit in phase 1
init_vector a(2,3,2) .. fit in phase 2

mfexp( a1 + a(2)*w + a(3)*w^2 )


Thanks guys. It runs well now with both methods you mentioned. I seems
like it is not good to have perfect data! I gonna try the random effect
module now...

Cheers,

Sylvain

2009/4/9 Mark Maunder <mmau...@iattc.org>

I see that Dave found that it was related to having no error in the
simulated data. I guess if you fixed the standard deviation (use a
negative phase) it might work without adding error.

Mark

From: Sylvain Bonhommeau [mailto:sylvain.b...@gmail.com]
Sent: Thursday, April 09, 2009 8:39 AM
To: Mark Maunder
Subject: Re: [ADMB Users] "Hessian does not appear to be positive
definite"

Hi Mark,

Thanks for your answer.
I simulated the data with R (code below). There is no error added to
data for now but that is the aim of today. I specified a pin file and
the initial conditions are not (but close to) the "known" parameters.

Thanks,

Sylvain

Sylvain Bonhommeau

unread,
Apr 9, 2009, 11:44:22 AM4/9/09
to da...@otter-rsch.com, us...@admb-project.org
Hi,

Thanks for your answer. That's what I did or tried to do and that's the results for the dataset with known parameters.  Here are the .tpl, .pin and .dat files.

Thanks a bunch,

Sylvain

2009/4/9 dave fournier <ot...@otter-rsch.com>
kd.tpl
kd.pin
kd.dat

Sylvain Bonhommeau

unread,
Apr 23, 2009, 2:12:35 PM4/23/09
to da...@otter-rsch.com, us...@admb-project.org
Hi,

I want to test the ability of my procedure to estimate parameters so I simulated data with known parameters. As you mentioned previously, when there is no errors, it is not possible to get the standard deviation. However in the .par file, the value of parameters are almost perfect. The problem is when I add a small noise to my simulated data, i.e. multiplying by a random variable normally distributed with mean=1 and sd=0.01, the estimated parameters are not the one I chose...

I attach the .tpl, .pin and the two .dat file (without and with error) as well as the R script to simulate the data.

(this is a simplified version of the model I sent you before)


Thanks for your help,

Sylvain



2009/4/9 dave fournier <ot...@otter-rsch.com>
I think we should keep this on the list as these are issures whihc may
kd_model_test.tpl
kd_model_test.pin
kd_model_test_without_error.dat
kd_model_test_with_error.dat
data_generator.R

dave fournier

unread,
Apr 23, 2009, 4:10:56 PM4/23/09
to Sylvain Bonhommeau, us...@admb-project.org
Sylvain Bonhommeau wrote:

What you are coming up against here are the limits of estimability.
If you remove the std dev from the model and just fit the least squares
and then
calculate the mean-squared residuals you will get 0.000821022 so the model
fits the data almost perfectly. the fact that they are not the values
you used
to generate the data simply reflects the fact that there is not a lot of
information in the data
i.e. thee are different combinations of the model parameters which fit
the model (almost) equally well.

Dave

> Hi,
>
> I want to test the ability of my procedure to estimate parameters so I
> simulated data with known parameters. As you mentioned previously,
> when there is no errors, it is not possible to get the standard
> deviation. However in the .par file, the value of parameters are
> almost perfect. The problem is when I add a small noise to my
> simulated data, i.e. multiplying by a random variable normally
> distributed with mean=1 and sd=0.01, the estimated parameters are not
> the one I chose...
>
> I attach the .tpl, .pin and the two .dat file (without and with error)
> as well as the R script to simulate the data.
>
> (this is a simplified version of the model I sent you before)
>
>
> Thanks for your help,
>
> Sylvain
>
>
>
> 2009/4/9 dave fournier <ot...@otter-rsch.com

> <mailto:ot...@otter-rsch.com>>

> 2009/4/9 Mark Maunder <mmau...@iattc.org <mailto:mmau...@iattc.org>>

> Us...@admb-project.org <mailto:Us...@admb-project.org>
> http://lists.admb-project.org/mailman/listinfo/users

Reply all
Reply to author
Forward
0 new messages