Improve VectorSpace concept

[sdrave]

This pull request contains two major changes to the VectorSpace concept:

VectorSpaces are now the factories for new VectorArrays. Instead of instantiating the array directly, one now writes space.make_array(obj_to_be_warpped_as_vector_array). In particular, to implement a new VectorArray own now also subclasses VectorSpaceInterface. This new approach has the following benefits:

• VectorSpaces replace the subtype concept which was hard to understand by newcomers. In contrast, the idea of a VectorSpace which can produce vectors contained in itself should be very natural.
• VectorSpaces can have meaningful attributes (like dim) instead of a generic subtype attribute. Custom implementations of __repr__ are easily possible.
• By subclassing an existing VectorSpace, one can easily add additional data that is relavant in a given mathematical situation (e.g. the subdomain number for localized model reduction), and this data is then also attached (via the space) to all arrays created from the space.

Second, every VectorSpace has now a string id attribute that can be used to distinguish mathematically different vector spaces. The currently used default id is 'STATE' for all function spaces and 'SCALARS' for elements of R^n. In particular, the range of functionals and the source of vectors are tagged as 'SCALARS'. The same would apply to multidimensional input or output operators. This in particular allows to detect the 'kind' (#275, #310) of an Operator without any additional effort. As a consequence, the special functionals and vector_operators dicts of the Discretizations can go away. The same applies to the ss_operators, is_operators, so_operators and io_operators dicts in @pmli's code.

This pull request is still work in progress. In particular, i haven't updated the documentation yet. Before I proceed with the polish, @andreasbuhr, @pmli, @ftalbrecht, please let me know what you think!

[sdrave]

In particular, you may want to take a look at pymor.reductors.basic.

[pmli]

I like this idea very much! I thought at some point that it would make sense to have input, state, and output spaces, and having the id attribute makes this possible!

BTW, this reminds me that input/output space reduction (i.e. reducing the number of inputs and/or outputs) is also something people do.

[andreasbuhr]

Hi,

I do like this change.

A suggestion:
I think that the "scalars" class makes the code hard to read. From the context, the user has to guess that "scalars" is a space definition. And did you think about complex and real scalars? I suggest renaming "scalars" to "RealScalarsSpace".

[sdrave]

Ok, so after some more polishing, this should now be ready to merge. Since this is a pretty deep change, we should still expect a few bug popping up ..

@andreasbuhr, after discussion with @ftalbrecht, I decided to keep scalars as it is used in a lot of places and much quicker to write. I believe it should quickly be familiar to most users. Regarding complex dtypes, I would propose to add an additional attribute to NumpyVectorSpace (and scalars), to set a non-standard dtype. However, we should discuss this further in #259.

[sdrave]

Hi @pmli, you may also want to take a look at the sys-mor-with-new-vectorspaces branch, where I have updated your code to work with the new VectorSpace concept (first commit).

I have then replaced the different ??_operators dicts by a single dict, using the source and range ids of the operators to identify their respective type. I used the SCALARS id to identify input and output spaces, which would not work if there were also state to input operators. For this, one should use 'INPUT' and 'OUTPUT' ids. However, using 'SCALARS', the code might be more generic ..

Not having different operators dicts then allowed me to use DiscretizationBase as a base class for the iosys discretizations, which made their implementation much more simple (I hope ...). In the last commit I made a minor change which seem to make more mathematical sense to me. However, I am not an expert ..

Note that I only have tested if the heat.py and delay.py demos still work, so there might be a few bugs around.

Please let me know if you like the changes overall or if you feel that the new VectorSpace concept should undergo further changes to better fit your code. If you agree with the overall direction I will merge this PR and open a new PR for the sys-mor-with-new-vectorspaces branch.

[pmli]

Hi @sdrave. I like the changes you made. I don't see anything right now that should be changed further in VectorSpaces.

Although, I must admit I'm not a big fan of using 'SCALARS' id for what could be higher-dimensional VectorSpaces, e.g. when there are multiple inputs and/or outputs. But, I'm not sure if by "scalars" you mean "a vector of scalars" or "elements of a field (or one-dimensional vector space)". In the first post here, you said that 'SCALARS' id is "for elements of R^n" and also for "the range of functionals and the source of vectors" 😕

State to input (or output to input) operators could appear in closed loop systems, and there are MOR methods for them, but I haven't looked into them very much. Maybe we can postpone this for later, I think you can open a new PR for the sys-mor-with-new-vectorspaces branch.

[sdrave]

Ok, so 'SCALARS' generally seems to be confusing.

To make my original intention clearer: for me, as a PDE guy, there is usually a big difference between state space (a discrete function space) and (possibly) multi-dimensional input/output spaces, which are typically something like R^n, where n is a relatively small number. The distinction whether the input space is one or three-dimensional usually is not as important as whether it is an Euclidean space or a function space. So I would have been fine with using 'SCALARS' for all those spaces (every field is a one-dimensional vector space over itself).

Having considered this a bit more, I would say that this special id is maybe not as needed as I first thought. So, how about defaulting the VectorSpace id to None instead of 'STATE', and only assign a special id, e.g. 'STATE', where useful?

Only one drawback, which was actually one of the motivations for introducing 'SCALARS': as a next step, I would like to add an as_array method to OperatorInterface which generalizes as_vector to multi-dimensional source or range spaces (the goal is to replace the op.apply(op.source.from_data(np.eye(op.source.dim))) idiom in @pmli's code). However, without 'SCALARS', there is no way to decided if the array should be from the source or from the range space. So there would have to be as_source_array and as_range_array methods. I guess, this is fine, as well.

[pmli]

Ok, the good thing is that in systems and control theory there is also a distinction between state, input, and output. So, the 'STATE' id seems to make sense in general (if we ignore the behavioral approach, which handles things differently, and state space formulations anyway seem more prevalent in practice).

I would definitely like to have a more general version of as_vector in the OperatorInterface. At first, I thought you could implement as_array by checking which space id is different from 'STATE', but then I realized that nothing would work for input to output operators (the D matrix in LTISystem), so I guess it would be good to have both as_source_array and as_range_array (or add an argument to as_array).

Maybe None as the default id would work fine. For interconnected systems (different systems connected through inputs and outputs), I could imagine having 'STATE1', 'STATE2', 'INPUT1', 'INPUT2', etc.

[sdrave]

As discussed, I have removed the special 'SCALARS' id and defaulted the id to None. Thus, there is no longer a notion of 'low-dimensional vector space which is not a function space' in pyMOR.

To check if a space is a one-dimensional NumpyVectorSpace, I have added the is_scalar attribute to VectorArrayInterface. This is used in pyMOR's discretizations to check if the given operators are of the right form. Note that we have space.is_scalar == (isinstance(space, NumpyVectorSpace) and space.dim == 1), whereas space.is_scalar != (space == NumpyVectorSpace(1)) in general, since space could have an id such as 'INPUT'.

I have implemented the as_source_array and as_range_array methods discussed above, which should be used in favor of as_vector, which we might want to deprecate. For now, I have fixed as_vector by adding an optional space parameter which needs to be specified for the corner case that op.source.is_scalar and op.range.is_scalar and op.source.id != op.range.id (this actually happens regularly when the RB space is one-dimensional).

To support, say, multi-input to state operators in estimate_image, I have disallowed functionals to be specified in the vectors argument (otherwise, detecting the space in which the image lies would have become quite nasty). Thus, if you have a functional, you now have to specify the transposed operator. To facilitate this, I have added the T property to OperatorInterface which gives you the transposed operator.

Note that I am speaking of 'transposed' and not 'adjoint' operator: I currently have the impression that apply_adjoint is over-engineered and should loose its product parameters. But then we have the transposed operator and not the adjoint operator in general, so we might as well rename the method to apply_transpose. I will open a separate issue for this.

[sdrave]

Thinking of the various corner cases I had to consider, I still feel that having something like 'SCALARS' would be nice. The straightforward solution to this would be to subclass NumpyVectorSpace. However, I am still wondering what the right semantics for such spaces would be in pyMOR, and I guess from @pmli's POV this might be even less clear. A definition like 'a low-dimensional VectorSpace which has a special meaning in some methods' is not really that convincing. Any more input on this, @pmli, @ftalbrecht, @andreasbuhr?

Implementing something such special spaces later on should be very easy. So I would vote for merging the code in its current state.

[pmli]

I have difficulty finding places where something like 'SCALARS' would be nice. @sdrave Could you point me to some examples?

I'm also not sure what you mean by 'function space', 'low-dimensional function space', and 'low-dimensional Euclidean space'. Is the first something like L^2(spatial domain), the second some subspace of the first (e.g. a finite element discretization), and the third the equivalent R^n to the second? Do each of these spaces have a representation in pyMOR?

[sdrave]

Sorry for the delay! @pmli, I have thought some more about this and still have no clear concept of what a 'scalar' should be. There are two places where it would be nice to have something like that:

• As I explained above, as_vector is bound to fail if source and range are one-dimensional but different (e.g. have different id). With something like 'SCALARS', it would be reasonable to assume that as_vector only works in cases where either source or range is a scalar VectorSpace. Making this assumption, as_vector would always work without specifying an additional space. Similarly, one could have an as_array method instead of as_source_array and as_range_array. Of course, one could also want to call as_array for Operators which map between two scalar spaces (e.g. input to output), in which case as_array again would be unable to decide whether to return the matrix or its transpose.

• For estimate_image, it is reasonable to assume that one wants to determine the space of possible vector representations of a single given parametric functional or vector-like operator. Without specifying a VectorSpace in which the image is to be determined, one is again left with the issue to determine whether the given Operator represents a functional or a vector (resp. an output or input Operator). I have solved this by only accepting vector-like operators, forcing the user to specify the transpose of functionals. With the concept of scalars, it would be possible for the algorithm to detect if the Operator should be seen as a functional or as a vector. Still, the issue of intput to output operators would remain.

Overall, however, I am fine with the current situation, so from my POV this is ready to be merged.

[pmli]

I also ok with merging. Maybe I'll understand things better after using it.