ASA with IDAS

[jbazil]

Hello,

I have a question regarding the use of IDAS to evaluate ~40 parameter sensitivities to several model outputs (consisting of both steady state and transient points).  Some of the model outputs are simply the state variables; however, some are functions of the state variables and parameters. There are a ~100 state variables with ~10 of those being algebraic variables.  I am using SundialsTB since the original code was written in MATLAB.  I do not have an analytical expression for the Jacobian.  I am thinking of playing with an AD algorithm to evaluate the Jacobian, but that's just a side note.   I have successfully run IDAS without the sensitivity modules active on the system of DAEs.  What is the best approach for the sensitivity analysis?

The forward sensitivity analysis is straight forward, but the size of the problem makes this approach unattractive.  I'm thinking the adjoint sensitivity analysis is the best way to go, but I am not sure exactly how to implement it.  Do I need to analytically solve for the adjoint equations?  Or can Sundials do that if given f(t,x,x',p)?  Also, I am not quite clear how to evaluate parameter sensitivities with respect to a model output using the adjoint method.  Do I need to solve for the analytical expressions for dg/dp to go into some residual function?  Or can Sundials generate these via finite difference approximations if I give g(x,p)?  In order evaluate g(t,x,p), I'd have to run several ASAs consisting of evaluating g(t1,x,p), g(t2,x,p), etc.?  To sum it all up, what are the minimum inputs needed for Sundials to run an ASA on g(x,p)?

I could not find an example problem that evaluated the parameter sensitivities for a vector valued objective function.  After some searching, it seems that only scalar valued objective functions are allowed, and I'd have to rerun the backwards integration scheme with a new scalar valued objective function.  Is this correct?  Any help would be greatly appreciated.  Thanks.

-jason

[Alan Hindmarsh]

Hello Jason,

I am not the expert here on Sensitivity Analysis (Radu Serban is).
But here is what I believe are the answers to your questions:

There is nothing preventing the objective functional integrand g (and
hence its integral G) from being a vector.  Just read Section 2.6 in
the User Guide with that vector interpretation.  The size of the
adjoint vector lambda is (size of y)*(size of g).

As things stand, you do have to supply the backward residual functions
(for the adjoint equations) analytically, and same for dg/dp and dF/dp
in the backward quadratures.

In principle, you should be able to get all desired results with a single
forward run and a single backward run.  But IDAS allows for solving multiple
backward problems for the same forward problem, if that is preferable
(e.g. to keep storage to a minimum).

If any of this is wrong, I trust Radu will correct me.

-Alan H


On Thu, 9 Jul 2009, jbazil wrote:

>
> Hello,
>
> I have a question regarding the use of IDAS to evaluate ~40 parameter
> sensitivities to several model outputs (consisting of both steady state and
> transient points).  Some of the model outputs are simply the state
> variables; however, some are functions of the state variables and
> parameters. There are a ~100 state variables with ~10 of those being
> algebraic variables.  I am using SundialsTB since the original code was
> written in MATLAB.  I do not have an analytical expression for the Jacobian.
> I am thinking of playing with an AD algorithm to evaluate the Jacobian, but
> that's just a side note.   I have successfully run IDAS without the
> sensitivity modules active on the system of DAEs.  What is the best approach
> for the sensitivity analysis?
>
> The forward sensitivity analysis is straight forward, but the size of the
> problem makes this approach unattractive.  I'm thinking the adjoint
> sensitivity analysis is the best way to go, but I am not sure exactly how to
> implement it.  Do I need to analytically solve for the adjoint equations?
> Or can Sundials do that if given f(t,x,x',p)?  Also, I am not quite clear
> how to evaluate parameter sensitivities with respect to a model output using
> the adjoint method.  Do I need to solve for the analytical expressions for
> dg/dp to go into some residual function?  Or can Sundials generate these via
> finite difference approximations if I give g(x,p)?  In order evaluate
> g(t,x,p), I'd have to run several ASAs consisting of evaluating g(t1,x,p),
> g(t2,x,p), etc.?  To sum it all up, what are the minimum inputs needed for
> Sundials to run an ASA on g(x,p)?
>
> I could not find an example problem that evaluated the parameter
> sensitivities for a vector valued objective function.  After some searching,
> it seems that only scalar valued objective functions are allowed, and I'd
> have to rerun the backwards integration scheme with a new scalar valued
> objective function.  Is this correct?  Any help would be greatly
> appreciated.  Thanks.
>
> -jason
> --
> View this message in context: http://*www.*nabble.com/ASA-with-IDAS-tp24412438p24412438.html
> Sent from the Sundials - Users mailing list archive at Nabble.com.
>
>
>

[Radu Serban-3]

Jason,

Sorry for the late reply... Alan's comments are all correct; I only have a couple of clarifications:

1. In principle, one adjoint problem solution (involving a forward pass and a backward one) is required to find the gradient with respect to the problem parameters for each output functional of interest. These outputs can be either in integral form G(p)=\int_t0^t1 g(t,x,p) dt or simply the "point-wise" g(t1,x(t1),p). If I understand it correctly, you are interested in the later type, but you have several different functionals g_i and you want their sensitivities at several different end-times t1_j. As such, you will need the solution of a (different) adjoint problem for each combination {i,j}.

While IDAS can take advantage of the fact that all these adjoint problems depend on the same forward problem (and, as Alan mentioned, solve all of them using a single backward pass), depending on how many functionals g you have and how many times t1 you are interested in, the cost of calculating these sensitivities with ASA may be in fact as large (or worse) than that of using FSA.

2. You do have to provide the adjoint problem corresponding to each output functional. This means you must be able to evaluate all partial derivatives, both for F and g, with respect to both x and p, either analytically or using (reverse) AD. Note that dg/dp (for some given functional g) is a vector of partial derivatives and should be, in principle, straight-forward to obtain. The terms in the adjoint problems that may be more difficult to obtain are in fact those

--Radu


Alan Hindmarsh wrote:

Hello Jason,

I am not the expert here on Sensitivity Analysis (Radu Serban is).
But here is what I believe are the answers to your questions:

There is nothing preventing the objective functional integrand g (and
hence its integral G) from being a vector.  Just read Section 2.6 in
the User Guide with that vector interpretation.  The size of the
adjoint vector lambda is (size of y)*(size of g).

As things stand, you do have to supply the backward residual functions
(for the adjoint equations) analytically, and same for dg/dp and dF/dp
in the backward quadratures.

In principle, you should be able to get all desired results with a single
forward run and a single backward run.  But IDAS allows for solving multiple
backward problems for the same forward problem, if that is preferable
(e.g. to keep storage to a minimum).

If any of this is wrong, I trust Radu will correct me.

-Alan H


On Thu, 9 Jul 2009, jbazil wrote:

  

Hello,

I have a question regarding the use of IDAS to evaluate ~40 parameter
sensitivities to several model outputs (consisting of both steady state and
transient points).  Some of the model outputs are simply the state
variables; however, some are functions of the state variables and
parameters. There are a ~100 state variables with ~10 of those being
algebraic variables.  I am using SundialsTB since the original code was
written in MATLAB.  I do not have an analytical expression for the Jacobian.
I am thinking of playing with an AD algorithm to evaluate the Jacobian, but
that's just a side note.   I have successfully run IDAS without the
sensitivity modules active on the system of DAEs.  What is the best approach
for the sensitivity analysis?

The forward sensitivity analysis is straight forward, but the size of the
problem makes this approach unattractive.  I'm thinking the adjoint
sensitivity analysis is the best way to go, but I am not sure exactly how to
implement it.  Do I need to analytically solve for the adjoint equations?
Or can Sundials do that if given f(t,x,x',p)?  Also, I am not quite clear
how to evaluate parameter sensitivities with respect to a model output using
the adjoint method.  Do I need to solve for the analytical expressions for
dg/dp to go into some residual function?  Or can Sundials generate these via
finite difference approximations if I give g(x,p)?  In order evaluate
g(t,x,p), I'd have to run several ASAs consisting of evaluating g(t1,x,p),
g(t2,x,p), etc.?  To sum it all up, what are the minimum inputs needed for
Sundials to run an ASA on g(x,p)?

I could not find an example problem that evaluated the parameter
sensitivities for a vector valued objective function.  After some searching,
it seems that only scalar valued objective functions are allowed, and I'd
have to rerun the backwards integration scheme with a new scalar valued
objective function.  Is this correct?  Any help would be greatly
appreciated.  Thanks.

-jason
--
View this message in context: http://**www.**nabble.com/ASA-with-IDAS-tp24412438p24412438.html
Sent from the Sundials - Users mailing list archive at Nabble.com.



    

[jbazil]

Radu,

Thanks for your reply.  You did clear up some questions I had; however, I have
an additional, more basic question on the FSA in general.

I apologize if this is wasting your time, but my question is about integrating
all the Ns(1+Np) differential equations simultaneously for the FSA versus simply
applying a small perturbation to each parameter and integrating the 1+Np (1+2Np
if the centered difference quotient is used) differential equations.  Is it
simply error control that necessitates Ns(1+Np) equations to be solved
simultaneously?  

In other words, why can't you compare the unperturbed trajectory versus the
perturbed trajectory for each parameter?  I can see that this would necessitate
the computation of Np (2Np if using centered DQs) different Jacobians, while the
same Jacobian is used to integrate all the Ns(1+Np) equations.  Is this the
extra computations that is undesirable and make the FSA more computationally
efficient?

Perhaps I am fundamentally misunderstanding the meaning of the sensitivity
equations.  Any clarification on this matter would be greatly appreciated.  Thanks.

-jason

--
Jason Bazil, PhD Candidate
Weldon School of Biomedical Engineering
Purdue University


Quoting Radu Serban <radu@...>:

> Jason,
>
> Sorry for the late reply... Alan's comments are all correct; I only have
> a couple of clarifications:
>
> 1. In principle, one adjoint problem solution (involving a forward pass
> and a backward one) is required to find the gradient with respect to the
> problem parameters for each output functional of interest. These outputs
> can be either in integral form G(p)=\int_t0^t1 g(t,x,p) dt or simply the
> "point-wise" g(t1,x(t1),p). If I understand it correctly, you are
> interested in the later type, but you have several different functionals
> g_i and you want their sensitivities at several different end-times
> t1_j. As such, you will need the solution of a (different) adjoint
> problem for each combination {i,j}.
>
> While IDAS can take advantage of the fact that all these adjoint
> problems depend on the same forward problem (and, as Alan mentioned,
> solve all of them using a single backward pass), depending on how many
> functionals g you have and how many times t1 you are interested in, the
> cost of calculating these sensitivities with ASA may be in fact as large
> (or worse) than that of using FSA.
>
> 2. You do have to provide the adjoint problem corresponding to each
> output functional. This means you must be able to evaluate all partial
> derivatives, both for F and g, with respect to both x and p, either
> analytically or using (reverse) AD. Note that dg/dp (for some given
> functional g) is a vector of partial derivatives and should be, in
> principle, straight-forward to obtain. The terms in the adjoint problems
> that may be more difficult to obtain are in fact those
>
> --Radu
>
>
> Alan Hindmarsh wrote:
> > Hello Jason,
> >
> > I am not the expert here on Sensitivity Analysis (Radu Serban is).
> > But here is what I believe are the answers to your questions:
> >
> > There is nothing preventing the objective functional integrand g (and
> > hence its integral G) from being a vector.  Just read Section 2.6 in
> > the User Guide with that vector interpretation.  The size of the
> > adjoint vector lambda is (size of y)*(size of g).
> >
> > As things stand, you do have to supply the backward residual functions
> > (for the adjoint equations) analytically, and same for dg/dp and dF/dp
> > in the backward quadratures.
> >
> > In principle, you should be able to get all desired results with a single
> > forward run and a single backward run.  But IDAS allows for solving
> multiple
> > backward problems for the same forward problem, if that is preferable
> > (e.g. to keep storage to a minimum).
> >
> > If any of this is wrong, I trust Radu will correct me.
> >
> > -Alan H
> >
> >
> > On Thu, 9 Jul 2009, jbazil wrote:
> >
> >  
> >> Hello,
> >>
> >> I have a question regarding the use of IDAS to evaluate ~40 parameter
> >> sensitivities to several model outputs (consisting of both steady state
> and
> >> transient points).  Some of the model outputs are simply the state
> >> variables; however, some are functions of the state variables and
> >> parameters. There are a ~100 state variables with ~10 of those being
> >> algebraic variables.  I am using SundialsTB since the original code was
> >> written in MATLAB.  I do not have an analytical expression for the
> Jacobian.
> >> I am thinking of playing with an AD algorithm to evaluate the Jacobian,
> but
> >> that's just a side note.   I have successfully run IDAS without the
> >> sensitivity modules active on the system of DAEs.  What is the best
> approach
> >> for the sensitivity analysis?
> >>
> >> The forward sensitivity analysis is straight forward, but the size of the
> >> problem makes this approach unattractive.  I'm thinking the adjoint
> >> sensitivity analysis is the best way to go, but I am not sure exactly how
> to
> >> implement it.  Do I need to analytically solve for the adjoint equations?
> >> Or can Sundials do that if given f(t,x,x',p)?  Also, I am not quite clear
> >> how to evaluate parameter sensitivities with respect to a model output
> using
> >> the adjoint method.  Do I need to solve for the analytical expressions
> for
> >> dg/dp to go into some residual function?  Or can Sundials generate these
> via
> >> finite difference approximations if I give g(x,p)?  In order evaluate
> >> g(t,x,p), I'd have to run several ASAs consisting of evaluating
> g(t1,x,p),
> >> g(t2,x,p), etc.?  To sum it all up, what are the minimum inputs needed
> for
> >> Sundials to run an ASA on g(x,p)?
> >>
> >> I could not find an example problem that evaluated the parameter
> >> sensitivities for a vector valued objective function.  After some
> searching,
> >> it seems that only scalar valued objective functions are allowed, and I'd
> >> have to rerun the backwards integration scheme with a new scalar valued
> >> objective function.  Is this correct?  Any help would be greatly
> >> appreciated.  Thanks.
> >>
> >> -jason
> >> --
> >> View this message in context:
> http://***www.***nabble.com/ASA-with-IDAS-tp24412438p24412438.html
> >> Sent from the Sundials - Users mailing list archive at Nabble.com.
> >>
> >>
> >>
> >>    
> >
> >  
>
>

[Radu Serban-3]

Jason Bazil wrote:

> Radu,
>
> Thanks for your reply.  You did clear up some questions I had; however, I have
> an additional, more basic question on the FSA in general.
>
> I apologize if this is wasting your time, but my question is about integrating
> all the Ns(1+Np) differential equations simultaneously for the FSA versus simply
> applying a small perturbation to each parameter and integrating the 1+Np (1+2Np
> if the centered difference quotient is used) differential equations.  Is it
> simply error control that necessitates Ns(1+Np) equations to be solved
> simultaneously?  
>
> In other words, why can't you compare the unperturbed trajectory versus the
> perturbed trajectory for each parameter?  I can see that this would necessitate
> the computation of Np (2Np if using centered DQs) different Jacobians, while the
> same Jacobian is used to integrate all the Ns(1+Np) equations.  Is this the
> extra computations that is undesirable and make the FSA more computationally
> efficient?
>
> Perhaps I am fundamentally misunderstanding the meaning of the sensitivity
> equations.  Any clarification on this matter would be greatly appreciated.  Thanks.
>
> -jason
>
>  

Jason,

I'm not quite sure I understand where you got the Ns(1+Np). In the
CVODES and IDAS user guides, we use Np to represent the total number of
problem parameters and Ns to represent the number of those problem
parameters with respect to which sensitivities are to be computed (hence
Ns<=Np). However, we say that the size of the combined system integrated
for FSA is N(1+Ns) or, if we will be computing sensitivities with
respect to all problem parameters, N(1+Np). Here N is the size of the
original problem. In other words, the cost of performing FSA is roughly
(1+Np) the cost of integrating the original system. So, from this
perspective alone, the computational cost is the same as if you were to
estimate the solution sensitivities using forward finite differences (in
which case, as you mention, you'd need to integrate the original
unperturbed system and then Np additional perturbed systems, each one of
these of size N).

Calculating 1st order solution sensitivities using continuous FSA vs.
estimating them with FD is more accurate, more efficient, and does not
require coming up with an appropriate perturbation factor.

Hope that this helps,
--Radu