ADMB sum-to-zero constraint and penalty for deviations

[Bai Li - NOAA Affiliate]

Hi all, 

I am implementing a recruitment module in Fisheries Integrated Modeling System with Ian and Kyle (cc'd) and trying to add a sum to zero constraint to recruitment deviations. We thought it would be great to learn from ADMB implementation, but we have a few questions after looking at the source code online. 

1) Does ADMB sum deviations to zero internally by subtracting the mean(vector of deviations) from the vector of deviations? Where is the source code for the implementation? Is it related to param_init_boudded_dev_vector?

2) Does ADMB implement an internal penalty associated with sum-to-zero constraint to make sure the sum-to-zero vector of deviations maps back to a unique set of deviations in the estimation space? Where is the source code for this penalty and what is the order of steps? I am getting a little bit lost by looking at ::set_value(*this,x,ii,minb,maxb,pen). Does the penalty show up in the NLL later?

It would be great if ADMB developers can help us better understand the components and order of ADMB sum to zero constraint implementation? Thank you!

Best,

Bai 

--

Bai Li, Ph.D.

National Stock Assessment Modeling Team, Contractor with ECS in support of

NOAA Fisheries Office of Science & Technology | U.S. Department of Commerce

[Jim Ianelli - NOAA Federal]

great questions. One of the topics that Matthew was supposed to tackle back in the day (documenting some of the technical aspects of ADMB that were missing) but...

comments below

On Mon, Jun 13, 2022 at 11:03 AM 'Bai Li - NOAA Affiliate' via ADMB Users <us...@admb-project.org> wrote:

Hi all, 

I am implementing a recruitment module in Fisheries Integrated Modeling System with Ian and Kyle (cc'd) and trying to add a sum to zero constraint to recruitment deviations. We thought it would be great to learn from ADMB implementation, but we have a few questions after looking at the source code online. 

1) Does ADMB sum deviations to zero internally by subtracting the mean(vector of deviations) from the vector of deviations? Where is the source code for the implementation? Is it related to param_init_boudded_dev_vector?

I think it is but am unsure. model23.cpp in folder nh99 seems to have some code that may be related (copied the contents below).

2) Does ADMB implement an internal penalty associated with sum-to-zero constraint to make sure the sum-to-zero vector of deviations maps back to a unique set of deviations in the estimation space? Where is the source code for this penalty and what is the order of steps? I am getting a little bit lost by looking at ::set_value(*this,x,ii,minb,maxb,pen). Does the penalty show up in the NLL later?

I think this occurs as well but unsure how it exactly happens in the code. I remember that it was something that implementing it the penalty basically goes to zero and the condition is met (sum to zero). It could be that same code listed below? 

It would be great if ADMB developers can help us better understand the components and order of ADMB sum to zero constraint implementation? Thank you!

Best,

Bai 

--

Bai Li, Ph.D.

National Stock Assessment Modeling Team, Contractor with ECS in support of

NOAA Fisheries Office of Science & Technology | U.S. Department of Commerce

--
To unsubscribe from this group and stop receiving emails from it, send an email to users+un...@admb-project.org.

--

James Ianelli
NMFS/NOAA Building 4
7600 Sand Pt Way NE
Seattle WA 98115

 206 526 6510

Visit the ADMB project http://admb-project.org/

Also see Alaska groundfish stock assessments

[Dave Fournier]

This is an example of where two features combine to become a bug.

The first feature is a dev vector.  This is simply a vector whose components should sum to  zero.  So how does it get these components?   Well during the fitting phase

the function minimizer ("solver") is picking a big vector of parameters  x to test

by evaluating the objective function on  x.  So how to do this.  Well first the vector x

needs to be reexpressed in terms of the independent model parameters

                                        par1

                         x    --->   par2           

                                        ...

                                        par n

Now note the the par i can be numbers, vectors, or higher dimension arrays.Tthey are declared in the   PARAMETER_SECTION

and have the prefix init_    in their declaration. from a math point of view they are

the independent variables.   They are all derived from a base class which puts them on a list.  One can run through this list and do certain things to all of them.   the

things one wants to do are accomplished by invoking the desired virtual function.

Examples are set_value which goes from the x to the independent variable and

set_value_inv which goes in the other way.  there is also a function to count up

the components of vectors etc.

So for a dev vector the set value function is intended to ensure that the components

sum to zero.  This is done by subtracting off the mean of the x's (this is not all the x's but only the ones that go to this parameter vector) ie. 

               x -->  x-mean(x)  ---> dev_vector

Now the dev vector components will sum to zero independent of whatever the x sum to.  But the sum of the x's may get large and in any event the Hessian will be degenerate.  So a penalty  of the form penwt*sum(x)*sum(x) is added to the function being minimized.  In general one starts with a sum(x)=0 and it stays small.

The penalty is passed back from the set_value function. 

Now we add boundedness to the dev function.  What can happen is that some of the

components are (almost) tight on the bounds and the sum(x) is not quite zero.

When starting a new phase the sum(x) is subtacted from the x's so that the sum zero condiiton is satisfied exactly.  However now the constraints may be violated.

This may cause the gradient to change a lot. Hence two features dev and boundedness combine to form a bug.

[Richard Methot - NOAA Federal]

Here's a relevant SS3 issue discussion:

https://github.com/nmfs-stock-synthesis/stock-synthesis/issues/37

and somewhere in ADMB archives will be discussion of a shortcoming in its implementation.  Sum to zero does not get used in MCMC, so a draw that gets included in MCMC produces weird result when used in MCEVAL.

[William Stockhausen - NOAA Federal]

I suspect the following has been considered many times before: Taking a page from the mechanics (physics) literature, I think it's possible to express the effects of the sum-to-zero constraints in terms of Lagrangian multipliers (one for each constraint), with a corresponding quadratic term for each in the constraint term subtracted from original likelihood. That is, for each constraint the corresponding multiplier times the quadratic term (just the constraint squared) is added to the objective function. Each Lagrangian multiplier is then regarded as an additional parameter, and the minimization of the objective function is done with respect to the augmented parameter set (i.e., original + Lagrangian multipliers), where the gradient of the augmented objective function is zero wrto each of the original parameters and each of the Lagrangian multipliers. Mathematically, at least, the constraints are satisfied identically at the solution to the augmented minimization problem. Don't know what sort of problems one would run into from a numerical standpoint. But this wouldn't get you around the boundedness problem nor the MCMC issues.

William T. Stockhausen

Research Fishery Biologist, Alaska Fisheries Science Center

NOAA Fisheries | U.S. Department of Commerce

Office: (206) 526-4241

www.fisheries.noaa.gov

[Cole Monnahan - NOAA Federal]

The Stan team as considered this problem from a Bayesian perspective (i.e., amenable to specifying meaningful priors) too: https://mc-stan.org/docs/2_28/stan-users-guide/parameterizing-centered-vectors.html  -Cole

[Bai Li - NOAA Affiliate]

Thank you all for your great explanations and suggestions! Based on the discussion, for FIMS, perhaps we could implement ADMB approach for now with unit and integration tests. So we can verify if FIMS milestone 1 product produces similar results compared to existing stock assessments. We also have created a new issue on FIMS repository and documented the other great suggestions for future comparison work. 

Dave, thanks again for explaining the ADMB approach. Has the penalty penwt*sum(x)*sum(x) been used outside of the function by any chance? If so, could you direct me to the right source code. Thank you!

Best,

Bai

On Tue, Jun 14, 2022 at 4:00 PM Richard Methot - NOAA Federal <richard...@noaa.gov> wrote:

Here is an example of the situation that Dave describes.  I can see this happening more clearly now with the improved console output in ADMB 13.0.

Warning: the following parameters had issues
  recdev1[2016] is within 0.1% of lower bound: -3 < -3.00005 < 3

In this case, I was experimenting with narrowing up the bounds on the recdev1 vector and when the model finished running, one of the devs has been rescaled to end up slightly outside the bound.

Rick

On Mon, Jun 13, 2022 at 2:54 PM Mark Maunder <mmau...@iattc.org> wrote:

Hi everyone,

 

You might want to consider if it is the right thing to do. I have not thought about it too much, but if your Rdevs don’t sum to zero and you have the standard penalty, then your model is probably misspecified (or using penalized likelihood causes bias). E.g. You could be using R0 incorrectly in creating the initial conditions. Also, a better approach to using R0 in management quantities might be to use average R (adjust appropriately by the SR relationship) over the appropriate period.

 

Also note that the R0 profile diagnostic performs differently whether you use sum to zero or not.