Confirm consistent usage of getReferenceToJacobian

[jtencer]

If the Jacobian (decoder not ODE) is not constant, we need to confirm that the correct Jacobian is used in all cases. This would also enable a time-varying Jacobian.

[fnrizzi]

@jtencer can you please list here what are the rules in your view that one should check for making sure the "correct jacobian" is used?

[jtencer]

Currently, the user defined mapping interface includes getReferenceToJacobian() which is called whenever the Jacobian is queried. To support nonlinear mappings and time-varying Jacobians, we need to:

1. Confirm that there is no scenario when the Jacobian could need to be updated in between calls to getReferenceToJacobian
2. Include the ROM state as an argument to getReferenceToJacobian i.e. jacobian_type & getReferenceToJacobain(const rom_state_type & romState)
3. Confirm that getReferenceToJacobian is only called once per rom state update to avoid unnecessary updates.

Because 2 has not been completed, my hacky implementation currently updates the Jacobian as part of applyMapping which results in unnecessary updates.

[fnrizzi]

@jtencer @eparish1 @pjb236

Ok, let's go in order, since this is important:

• the decoder's Jacobian is taken by reference. so if you change it inside the custom decoder at any time, this will affect pressio too. This is dangerous but also helps avoiding copying data back and forth. So i will leave it like this: you are in charge of storing the decoder jacobian, and pressio only views it.

• I will fix pressio to make sure it calls getReferenceToJacobian only once to store the reference to make things more clear. This should not change anything as far as the numerics for the current code because calling getReference many times, simply restores the references which is ugly, so i will fix this. Since pressio only calls getReferenceToJacobian once, that method shoudl stay without arguments because the whole purpose of that method is to store the reference.

• so, now, my question is: for time-varying decoder's Jacobians, you want to recompute that Jacobian once for every rom state update, correct? If this is true, then the decoder class should have another method:

  // something like
  void updateJacobianIfNeeded(romState, time){}
  // or (we can find the best name)
  void updateJacobian(romState, time){}

  which pressio calls only once after every rom state update. So basically if you need to, you update the decoder's Jacobian once and that Jacobian is used for computing both the LSPG/WLS/Galerkin operators. Agree?

• by "rom state update" do you mean a full iteration, not the intermediate steps during the nonlinear solver. Right?

[jtencer]

In the regression test I have void updateJacobian(const rom_state_type & romState) const

My thought was that we would call that method from getReferenceToJacobian rather than having Pressio call it separately. This was the reason I wanted to avoid spurious calls to getReferenceToJacobian. If Pressio calls updateJacobian separately, then additional calls to getReferenceToJacobian aren't a big deal because there would be no computation involved in that op.

I think it's appropriate to only update the Jacobian once per full iteration. I think that if the Jacobian is changing significantly within an iteration, that either the step is too large or the manifold is too crinkly.

[fnrizzi]

@jtencer @eparish1 @pjb236

• what if I call the updateJacobian method before starting a new step? so basically, the order for LSPG would be at each time step: (a) update decoder Jacobian, (b) compute ROM residual, (c) compute ROM Jacobian, (d) solve least-squares problem.

• one problem wtih changing Jacobian: if the Jacobian changes at every time step, do we need to reconstruct the FOM states for the ODE stencil with the new jacobian or not? right now, at every time step, I only reconstruct the current and the previous FOM state. All other are just "replaced" by removing the oldest. But if we update the Jacobian at every time step, this might be a problem since some FOM states in the ODE stencil might be those reconstructed with previous mappings.

[jtencer]

I think we need to write out the equations we’re solving and figure out at what level we’re ok with linearization. That should be what governs the frequency we update the Jacobians.

[fnrizzi]

Yes, makes sense... would you be able, please, to start this? and then we can follow up here. This being said, I will not include this in the new upcoming release ok? :)

[jtencer]

I think we need to call computeJacobian each time we compute the ROM Jacobian. I don't think we need to recompute the old steps though.

[fnrizzi]

This relates to #107

I think we need to call computeJacobian each time we compute the ROM Jacobian. I don't think we need to recompute the old steps though.

Ok. So why don't we need to recompute old steps? can you post a pdf here with the math you did to figure this out?

[jtencer]

Because the appropriate decoder Jacobian for computing the old steps is the one that's evaluated at the old step. The time integrator shouldn't need to be aware of anything besides the history of states and the method for computing the next state given the prior states. What would be the justification for changing prior state information based on the current state?

[fnrizzi]

What would be the justification for changing prior state information based on the current state?

None! I agree with you! just wanted to make sure :) I will work on it for the new release

[fnrizzi]

@jtencer @eparish1 @pjb236 I am working on this and all is fine for explicit Galerkin.
But came across a doubt when working on implicit Galerkin with the discrete-time API (but I think the same happens with LSPG). Let's say we are a at a newton step k and let's call the decoder's jacobian J_d.

If we update the decoder's Jacobian ONLY when we compute the ROM Jacobian, then we have a discrepancy because at k the ROM residual computation will use the previous J_d while the computation of the ROM Jacobian will use the new J_d. So shouldn't we update the decoder's Jacobian once at the beginning of every solver's step?

[pjb236]

I don't see a way around that if we want to be mathematically correct. J_d should be computed whenever we update the state. However, updating J_d for every Newton step could be very costly.

[fnrizzi]

what if pressio passes you also the time step we are at? do you need the solver iteration too? that way you can decide whether to update or not. would that make sense?

  template<typename gen_coords_t>
  void updateJacobian(const gen_coords_t & genCoordinates) const
  {
  }

if you have a linear decoder, then of course the above is a noop but for a generic mapping it makes sense right?

[fnrizzi]

• solver owns when to update Jacobian

• updateJacobianFrequency updateJacobianEveryNSteps : for nonlinear solver Jacobian is J