Towards linear optimization (dual problem, sensitivity problem, output gradient)

[TiKeil]

This PR adds code for handling PDE-constrained parameter optimization. We only cover simple linear optimization in the sense that the output functional is linear. For that, we add code to compute the gradient of the output functional w.r.t. the parameter.

To compute the gradient, two approaches are feasible. One is to compute the sensitivities of the solutions and one is the adjoint approach, where a dual solution is used.

We also added a demoscript where we test all of the code by comparing the computational time between FOM and ROM implementation.

For more information on the topic, we also have a tutorial in #1205,

[sdrave]

Not complete sure about the interface myself, but some remarks:

• I would rename dual_model to dual. What else should the dual of a model be if not a model.
• I would ditch solve_dual completely. The user simply can write m.dual.solve()
• I don't like the name solution_sensitivity, but I'm not sure about a better name. How about solve_d_mu?
• I would say that the actual computation of solution_sensitivity should happen in Model.compute. I'm going to prepare a PR for the discussed changes tomorrow.
• output_gradient sounds better than solution_sensitivity. But maybe also output_d_mu?
• I would not add P and adjoint_approach to the signature in Model. I would assume this does not make much sense for arbitrary Models.
• One could consider base class implementations via finite differences.

[TiKeil]

Thanks, @sdrave for your very nice input.

Not complete sure about the interface myself, but some remarks:

* I would rename `dual_model` to `dual`. What else should the dual of a model be if not a model.

* I would ditch `solve_dual` completely. The user simply can write `m.dual.solve()`

I agree !

* I don't like the name `solution_sensitivity`, but I'm not sure about a better name. How about `solve_d_mu`?

I think this is a good name. I agree

* I would say that the actual computation of `solution_sensitivity` should happen in `Model.compute`. I'm going to prepare a PR for the discussed changes tomorrow.

Ah, even better !!

* `output_gradient` sounds better than `solution_sensitivity`. But maybe also `output_d_mu`?

agreed

* I would not add `P` and `adjoint_approach` to the signature in `Model`. I would assume this does not make much sense for arbitrary `Models`.

I think an adjoint approach is always possible for computing the gradient arbitrary outputs of arbitrary models but maybe that goes a little bit too far. So I agree to removing it there.

* One could consider base class implementations via finite differences.

That is a very nice idea.

[TiKeil]

I adjusted this code to the new _compute methodology. There are some things that I do not like yet.

1. I use the name output_d_mu for the gradient of the output
2. I use the name solution_d_mu for the partial derivative to only one parameter component.
3. I use the name solution_d_mus for a dict with all parameter components.

In my opinion the names of 1. and 2. do not really fit to each other, but so far I did not find a better name for these methods.

Apart from that, I think this PR is a good start for an optimization tutorial which I will prepare after this PR is merged.

[codecov]

Codecov Report

Merging #1110 (6abf6dd) into master (70af0f0) will increase coverage by 12.58%.
The diff coverage is 93.33%.

Impacted Files	Coverage Δ
src/pymortests/demos.py	88.48% <ø> (ø)
src/pymor/models/interface.py	79.69% <89.36%> (+13.40%)	⬆️
src/pymor/models/basic.py	89.89% <91.17%> (+12.97%)	⬆️
src/pymordemos/linear_optimization.py	96.20% <96.20%> (ø)
src/pymor/operators/block.py	84.72% <100.00%> (+0.36%)	⬆️
src/pymordemos/neural_networks.py	93.18% <0.00%> (ø)
src/pymordemos/parabolic_mor.py	96.26% <0.00%> (ø)
src/pymordemos/burgers.py	96.15% <0.00%> (ø)
src/pymordemos/thermalblock_adaptive.py	59.37% <0.00%> (ø)
... and 82 more

[sdrave]

How about extending compute to allow either passing a bool to solution_d_mu/output_d_mu or a tuple of parameter and index?

Maybe rename _compute_solution_d_mus to _compute_solution_d_mu to align it with the signature of compute and rename _compute_solution_d_mu to _compute_solution_d_mu to something like _compute_solution_d_mu_single_direction.

[sdrave]

Before merging this, we need to find a solution here. One simple option would be to let the user pass an optional dual_operator. However, I'm unsure if there are situations where also output_functional needs to be modified to be usable as a rhs.

[TiKeil]

I mean, it is not even clear that a dual model has the output_functional as right hand side. This probably is something that is usual for RB or Optimization, but in general one could also go for the dual model with the actual rhs as its right hand side.

We could thus also introduce an optional dual_rhs argument

[TiKeil]

Thanks @sdrave for the nice feedback. I have tried to address all your comments.

I am now using additional arguments dual_operator and dual_rhs for StationaryModel. This requires the user to add these arguments by hand once the discretizer has been called. While this might be a little bit annoying for the user, it is still better to ensure that the user indeed uses the correct rhs and operator for the dual model. Default quantities might be misleading.

PS: @renefritze Why did the enforce label check fail again?

[pymor-lab-hub-bridge]

Hey @TiKeil it looks like this PR touches our CI config (and you are not in my contributor safelist), therefore I am not starting a gitlab-ci build for it and it will not be mergeable.

[renefritze]

Hey @TiKeil it looks like this PR touches our CI config (and you are not in my contributor safelist), therefore I am not starting a gitlab-ci build for it and it will not be mergeable.

Yeah, shot myself in the foot there with the update, didn't I? 🙄

[renefritze]

Approving for the bridge bot to start the PR build again

[TiKeil]

@renefritze , something seems to be wrong with the wheels?

[renefritze]

@renefritze , something seems to be wrong with the wheels?

Yes

[TiKeil]

Thinking about it again, I think it is better to remain the question of "how to define a dual model?" open and open an issue for this instead. The reason is that there does not really exists a unique dual for a model because the right hand side might vary. Thus I hid the use of a dual model entirely in the computation of the gradient. For now, we do not need that somewhere else. For the future however, I would prefer to find a valid way for a incorporating a dual_model in the interface.

What is unique though is a dual_operator. If we do not find a better way of asserting that the operator is symmetric, we should still stick with dual_operator in the init of StationaryModel. @sdrave, If you agree with this way, we could merge this soon.

[sdrave]

@TiKeil, have you thought about my proposal of taking a discrete point of view and always using operator.H for the dual model?

[TiKeil]

@TiKeil, have you thought about my proposal of taking a discrete point of view and always using operator.H for the dual model?

Yea I did. operator.H can not be used immediately and results into a strange behavior with interesting solutions. So either one adapts .H or does not follow this.