On Thu, Aug 18, 2016 at 8:55 AM, Munson, Todd <tmun...@mcs.anl.gov> wrote:
> > Just to add some more thoughts to this argument. Currently TAO uses a > model where we have methods to get the gradient and Hessian of the > objective function. > > For quadratic problems we really want the Hessian matrix H and the > vector c so that the objective function is x'*H*x + c'*x. > > Because we do not have any specialized code for quadratic problems, > we end up using the gradient of the objective and the Hessian > matrix to back out the constant c. > > c = grad - Hx > > where grad = Hx + c. > > This relies on cancellation of Hx in order to get an accurate value > for the vector c, which can be a problem, especially if we start > far from a solution. > What about c = grad 0? Matt > You may have the same problem with linear ODE if you are relying > on a residual and jacobian evaluation pair to back out the > constant right-hand side vector. > > Its probably too complicated for users, but what one might like is > to have the user separate the objective function, for example, and > provide, c, H, and a nonlinear part f. Then I can assemble the > objective function internally as > > c'*x + x'*H*x + f(x) > > while retaining correct values for c and H when I need them. > > The user could continue as always by putting everything into f(x), > but then they would have to use a nonlinear solver. > > I'd do the same for constraints: > > b + A*x + F(x) > > This also allows the methods to do internal type checking to ensure > the problem passed into the methods has the correct signature. > E.g. for quadratic programs the nonlinear functions are NULL. > > Todd. > > > On Aug 16, 2016, at 7:03 PM, Barry Smith <bsm...@mcs.anl.gov> wrote: > > > > > >> On Aug 16, 2016, at 6:09 PM, Munson, Todd <tmun...@mcs.anl.gov> wrote: > >> > >> > >> Is there a reason to put flags into PETSc for linear or nonlinear > equations? In PETSc > >> there is a natural difference between linear (KSP) and nonlinear (SNES) > that is > >> reflected in the objects. > > > > The flags are for TS ODE solvers (not KSP or SNES); they can come into > play for two reasons > > > > 1) whether Jacobian entries can be changed inside the solver or not (for > nonlinear ODES they can usually be changed since the user provides new > values at each Jacobian update; for linear the user won't normally be > putting in the same values again so if the solver changes the values it has > to have some way to recover the values; often but not always the changes > are only on the diagonal). Note that this reason is JUST a question of > optimizing memory and code (do you make a copy of matrices or not?), not > correctness of results. > > > > 2) There may also be specialized ODE integrator only worked for a linear > problem. (I don't know if we have any of those implemented or not). This is > not just a question of optimizing the solver, it is a question of > correctness. If I try to solve a nonlinear problem with a linear only > solver I will get a wrong or meaningless answer. > > > > Similarly I think for TAO both reasons for wanting to know properties > of the problem exist. > >> > >> Any thoughts on splitting up the TAO methods in a similar way and > distinguish at > >> least between between quadratic and general nonlinear problems, rather > than > >> putting in flags. This would be consistent with the rest of PETSc. > >> > >> For example, instead of having a single "bound" directory, having > something > >> like a "bco" directory for generic nonlinear bound-constrained > optimization > >> and a "bqp" directory for bound constrained quadratic programs. > >> > >> I'd do the same for "unconstrained" into "uco" and "uqp". > >> > >> General constraints then become more cumbersome with too many > combinations. > >> Could probably get by with generic "lp", "qp", and "nlp" though. > > > > If you did this for the TAO user API then one could argue it should be > done for TS also. But I don't think it should be done for TS. > > > > Two reasons for not having a separate API for linear and nonlinear > problems is > > > > 1) I may want to use a "nonlinear" ODE integrator (that is an integrator > that works for nonlinear problems also) on my linear problem and I > shouldn't need to have ifs if my code to switch to use it. That is I don't > want ugly code like > > > > if (using nonlinear integrator) { > > TSSetIJacobian(..) > > } else { > > TSSetConstantMatrix(...) > > } > > Now I could also have my "nonlinear integrators" support the > TSSetConstantMatrix() API but that leads to more complex implementations > that are more work to maintain. Why not just have something like > > TSSetIJacobian(ts, J,TSConstantJacobianFunction()) when you are > doing a linear ODE? > > > > 2) I may spend 3 months writing my linear ODE problem using the linear > interface and then realize I have to add a nonlinear term, now I have to > change my app use the nonlinear interface of TS. This is why Matt argues we > should just publicly support the SNES api even for linear problems. > > > > Dealing with changing or not changing matrices has another issue we > have not dealt with properly in PETSc (or TAO). When the matrix entries are > set. For linear ODE problems we kind of assume the matrix entries are set > by the user before they ever call the TS code (and many users think this is > a natural way to proceed), so essentially they compute the J matrix up > front and then pass it to the TS which uses it. For nonlinear problems, of > course, we cannot set the matrix entries up front, they MUST be set in the > callback that TS calls to compute the matrix entries. But if the matrix > structure is provided internally by a DM then the user does not have direct > access to the matrix before calling TSSolve() so even for linear problems > the code logic more naturally wants the (constant) Jacobian to be computed > at the first timestep computation (and never again) rather than up front. > Using the Matt Knepley consistency/simplicity argument we should just > always use the callback approach and not even "support" the "linear only" > approach. Then comes the nasty linear "flag" to indicate 1) if the user > intends to provide the matrix only once and 2) to indicate the problem has > special (linear) structure that only certain solvers will work correctly on. > > > > I'm rambling a bit but I need to in order to make clear these kind of > design decisions are not straightforward. > > > > In terms of organizing the implementation code and the names (TaoType) > of the methods I think spitting them into categories makes good sense. But > never with your horrible abbreviations that no one but you would be able to > remember! We are not using Fortran IV with 6 char limits. > > > > Before doing anything let's try to generate a table that includes all > (most of?) the cases. > > > > The function being optimized is linear, quadratic, or general > (another case is separable). > > > > The equality constraints are linear?, quadratic? general? or > something else? > > > > The inequality constraints are simple bounds? linear? quadratic? > general? something else? > > > > Is there a fourth (or even more) axis of problem properties? > > > > Is this level of categories enough to organize the directories and > naming of TAO solvers? > > > > Barry > > > > > > > >> > >> The methods would then explicitly mention the type of optimization > problem > >> with, for example, TAO_UCO_LMVM, TAO_BCO_TRON, and TAO_BQP_GPCG. > >> > >> Thoughts? > >> > >> Todd. > >> > >> > >>> On Aug 16, 2016, at 4:32 PM, Oxberry, Geoffrey Malcolm < > oxber...@llnl.gov> wrote: > >>> > >>> Definitely flags for linear problems would be helpful for TAO. Once > there is an example up, I'd be happy to add that to the SQPTR pull request. > >>> > >>> ________________________________________ > >>> From: petsc-dev-boun...@mcs.anl.gov [petsc-dev-boun...@mcs.anl.gov] > on behalf of Barry Smith [bsm...@mcs.anl.gov] > >>> Sent: Monday, August 15, 2016 6:18 PM > >>> To: Munson, Todd > >>> Cc: petsc-dev > >>> Subject: Re: [petsc-dev] coding style > >>> > >>>> On Aug 15, 2016, at 8:08 PM, Munson, Todd <tmun...@mcs.anl.gov> > wrote: > >>>> > >>>> > >>>> Since we are only doing diagonal modifications to the matrix in > >>>> my case, should I simply create a matnest in my part of the > >>>> code and apply the diagonal modification in the matnest > >>>> without directly modifying any of the matrices? > >>>> That might be cleaner for the users. > >>> > >>> It will be dependent on the solver how much of the matrix you need > retain the parts you have modified. > >>> > >>> MatNest isn't the right thing here. I think you can do > MatGetDiagonal() to keep the valid diagonal entries then MatDiagonalSet() > to put back the valid entries. > >>> > >>> Barry > >>> > >>> > >>>> > >>>> Todd. > >>>> > >>>>> On Aug 15, 2016, at 7:55 PM, Barry Smith <bsm...@mcs.anl.gov> wrote: > >>>>> > >>>>> > >>>>>> On Aug 15, 2016, at 7:42 PM, Munson, Todd <tmun...@mcs.anl.gov> > wrote: > >>>>>> > >>>>>> > >>>>>> Got it. I will get all the code fragments in the developer doc that > >>>>>> are not meant to be deliberately wrong fixed tomorrow; at least > >>>>>> that code will be compliant... :) > >>>>>> > >>>>>> Then its onto fixing bqpip. > >>>>>> > >>>>>> Was there a final say on the Hessian/Jacobian matrix when they are > >>>>>> modified? Were we adding a "dirty" bit or was I changing my code > >>>>>> to not modify those matrices? > >>>>> > >>>>> Changing email to petsc-dev since this is a general PETSc/TAO issue. > >>>>> > >>>>> This is why we don't wait months to fix something :-) I have > completely forgot the context you are asking about. > >>>>> > >>>>> I think the general issue comes up when a problem is linear and the > user just wants to set and forget the matrices while when the problem is > nonlinear the user needs to reset values into the matrix anyways so it is > fine for you to change the matrix internally. Unless there is a way for the > user to explicitly indicate the system is linear your code cannot know if > it needs to make a copy of the matrix. > >>>>> > >>>>> In TS we have typedef enum {TS_LINEAR,TS_NONLINEAR} TSProblemType; > and TSSetProblemType() that can be used to indicate if the matrices can be > changed or not (though I don't think we use it in this way, we have an open > issue https://bitbucket.org/petsc/petsc/issues/135/memory- > optimization-for-ts > >>>>> > >>>>> I think you should start by having a flag that indicates if it is > safe to modify the matrices associated with the TAO object and if not safe > keep a backup copy. > >>>>> > >>>>> Barry > >>>>> > >>>>> > >>>>> > >>>>>> > >>>>>> Todd. > >>>>>> > >>>>>>> On Aug 15, 2016, at 6:09 PM, Barry Smith <bsm...@mcs.anl.gov> > wrote: > >>>>>>> > >>>>>>> > >>>>>>> Todd, > >>>>>>> > >>>>>>>> On Aug 15, 2016, at 10:26 AM, Satish Balay <ba...@mcs.anl.gov> > wrote: > >>>>>>>> > >>>>>>>> Todd, > >>>>>>>> > >>>>>>>> Forwarding this to Barry. > >>>>>>>> > >>>>>>>> Satish > >>>>>>>> > >>>>>>>> On Mon, 15 Aug 2016, Munson, Todd wrote: > >>>>>>>> > >>>>>>>>> > >>>>>>>>> I'm going through the developer documentation. There seems to be > >>>>>>>>> a few things missing from the coding style. > >>>>>>>>> > >>>>>>>>> For example, in function prototypes, is it: > >>>>>>>>> > >>>>>>>>> PetscErrorCode MyFunction(char *,const int *,int); > >>>>>>>>> > >>>>>>>>> or > >>>>>>>>> > >>>>>>>>> PetscErrorCode MyFunction(char*,const int*,int) > >>>>>>> > >>>>>>> The second one because the preference is to not have unnecessary > space. > >>>>>>> > >>>>>>>>> > >>>>>>>>> > >>>>>>>>> Same thing in the definitions: > >>>>>>>>> > >>>>>>>>> PetscErrorCode MyFunction(char *c,const int *i,int j) > >>>>>>>>> { > >>>>>>>>> } > >>>>>>>>> > >>>>>>>>> or > >>>>>>>>> > >>>>>>>>> PetscErrorCode MyFunction(char* c,const int* i,int j) > >>>>>>> > >>>>>>> The * should always be attached to the variable. > >>>>>>>>> { > >>>>>>>>> } > >>>>>>>>> > >>>>>>>>> The same questions apply with * replaced by []. > >>>>>>>>> > >>>>>>>>> > >>>>>>>>> Next is in the types for variables. Are pointers grouped > >>>>>>>>> separately or together with the type: > >>>>>>>>> > >>>>>>>>> int *i,j,*k; > >>>>>>>>> > >>>>>>>>> or > >>>>>>>>> > >>>>>>>>> int *i,*k; > >>>>>>>>> int j; > >>>>>>>>> > >>>>>>>>> If grouped together, should they be ordered with pointers first > >>>>>>>>> and the others later? > >>>>>>>>> > >>>>>>>>> int *i,*k,j; > >>>>>>> > >>>>>>> There is no preference on the above; the pointers do not have to > be on their own line. > >>>>>>>>> > >>>>>>>>> On the next topic, do you use > >>>>>>>>> > >>>>>>>>> char *blah; > >>>>>>>>> > >>>>>>>>> or > >>>>>>>>> > >>>>>>>>> char* blah; > >>>>>>>>> > >>>>>>>>> when there is only one of them? > >>>>>>> > >>>>>>> * belongs on the variable > >>>>>>> > >>>>>>>>> > >>>>>>>>> I have the same issues with the typedefs. Sometimes there is > >>>>>>>>> > >>>>>>>>> typedef struct _EH* EH; > >>>>>>>>> > >>>>>>>>> and other times > >>>>>>>>> > >>>>>>>>> typedef struct _EH *EH; > >>>>>>> > >>>>>>> Here the * belongs as in the first case > >>>>>>>>> > >>>>>>>>> > >>>>>>>>> My personal thought is to be consistent and use char * with a > space > >>>>>>>>> everywhere, but you tell me. I will then go and fix the > developer > >>>>>>>>> documentation to be consistent with the style guide as a test of > >>>>>>>>> my branch management. > >>>>>>>>> > >>>>>>>>> Todd. > >>>>>>>>> > >>>>>>>>> > >>>>>>>> > >>>>>>> > >>>>>> > >>>>> > >>>> > >>> > >> > > > > -- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener