Support complex data in NumpyVectorArray / NumpyMatrixOperator

[pmli]

The issue is that Numpy/Scipy doesn't allow a matrix to change type inplace, so a ComplexWarning: Casting complex values to real discards the imaginary part is raised in several places:

1. in NumpyMatrixOperator.assemble_lincomb in the line with

   matrix += (op._matrix * c)

   when the first coefficient is real and some later one is complex.
   Options here would be to use

   matrix = matrix + (op._matrix * c)

   or do a check at the beginning and initialize matrix appropriately (or always put complex coefficients first).

2. in pymor.operators.numpy._apply_inverse in the line with

   R[i] = np.linalg.solve(matrix, VV)

   The problem is that R is initialized with

   R = np.empty((len(V), matrix.shape[1]))

   An option here would be to add dtype=complex in np.empty if matrix or V is complex.

[ftalbrecht]

Just from a quick glance I would suggest to rather change the initialization than the operations, since there is often good reason to use += over = ... + .... Since we initialize new storage in both cases anyway (matrix and R) I would find it reasonable to create those as complex, if required.

[sdrave]

I tend to agree with @ftalbrecht that in your first example the best solution would be to have all matrices have complex dtype in the first place. Of course, this throws away memory and involves some performance penalty. @pmli, would doing so be a serious issue for you?

For the second example: I guess we need a more systematic approach for handling complex data. We could do so by adding the dtype to the subtype of NumpyVectorArray. However, this would again raise the issue if I am allowed to apply an operator with source subtype '(complex, dim)' to an array with source subtype '(float, dim)'. (See our discussion in #56.)

If we decide to add the dtype to the subtype (which sounds reasonable to me), there seem to be two ways to go:

1. Relax the notion of being contained in a VectorSpace by letting the array implementation decide if vectors of subtype A are contained in the space with subtype B. However, it is not clear to me what the exact semantics should be.
2. Add methods for casting arrays of given subtype into arrays of another subtype. Downside here would be that this would probably lead to unneeded copies of data.

[ftalbrecht]

Regarding the first example: I think the memory efficiency depends on numpys implementation of +=. It might be enough to leave all component matrices untouched and just create matrix as complex if one of the components is. Since we create matrix here anyway there should not be any wasted memory (if matrix += ... does the correct thing when the right hand side is real)...

[sdrave]

@ftalbrecht, I am not sure if I understand your comment. If operators[0]._matrix if stored as complex even though it is real, we need twice as much memory ...

[ftalbrecht]

Ah, sorry for not being clear enough! Consider the assemble_lincomb implementation of the NumpyMatrixOperator (link to docu):

All I propose is

• to detect, if any of the operators or coefficients is complex and
• then (and only then!) to instantiate matrix as a complex matrix
• to keep the rest (e.g., +=) as it was

If any of the operators is complex we need to create a complex matrix anyway, or the linear combination would not be correct any more; or am I wrong here? So we do not waste memory in that case. If none of the operators is complex we never instantiate a complex matrix, and do not waste any memory either. Or am I missing something here?

I do not propose to just store all matrices as complex (or convert them)...

[pmli]

I agree with @ftalbrecht, I would prefer to only use complex values when necessary.

@sdrave About dtype as subtype, multiplying a real matrix and a complex vector and vice versa both make sense, and both can occur in practice.

Also:
3. a number of operations in NumpyVectorArray don't work with complex numbers (scal, axpy, l2_norm...).

[sdrave]

Ok, so after some sleep, I agree that we probably should not overengineer things. If we get everything working with some additional checks, we should go that route. However, I am not sure if I am aware of all consequences. My proposal would therefore be to make some tweaks in the complex branch I have just created until @pmli's code is working well. If we are happy with the result, we merge to master.

[sdrave]

@pmli, pull requests are always appreciated. I do not have any code here working with complex numbers ...

[ftalbrecht]

Sounds reasonable, agreed!

[pmli]

Sounds good!

[pmli]

I've made some changes in the complex branch for NumpyMatrixOperator and NumpyVectorArray. Now I'm wondering about scalar products (dot and pairwise_dot). It makes some sense that a complex conjugation should appear somewhere. I would vote for a conjugation in the second position.

[sdrave]

I know that some libraries support both scalar products (with complex conjugation and without). As I do not see a lot of use for the variant without conjugation, I would say the conjugation should always be performed. And, yes of course, scalar products are antilinear in the second variable! 😀

[pmli]

Ok 😃
Unfortunately, not everybody agrees about which variable is supposed to be antilinear. Some have the first variable as antilinear, and I guess that then the Riesz representation theorem is "nicer".

While working with vector arrays, I noticed that it would be nice if they had real and imag attributes (like numpy arrays). Would you be ok with that?

[sdrave]

I think that having real and imag attributes would be nice. However, we should clearly define the semantics: Will accessing real on a float array return itself or a copy? (Don't know how NumPy behaves. Probable best to do it the same way.)

[pmli]

Looking at the examples in NumPy docs, it seems that real and imag always return a reference. Quick tests in the console confirm this. Therefore, changing real or imag changes the initial array.

[sdrave]

I have merged @pmli's changes in #138 which should do most of what is needed. We should leave this issue open until we have merged the complex branch into master. Before doing that we should take a look at the performance implications of the additional checks and optimize if needed.

[sdrave]

Regarding real and imag: I guess the safest bet would be to let real and imag always return a copy of the array. When #55 is implemented, we could use the same copy-on-write semantics for real and imag as well. @pmli, do you have use cases, where you would need a reference to the real or imaginary parts which can be modified?

[pmli]

@sdrave No, I don't have use cases when I need references (I never change real or imaginary part), and I would prefer to work with copies.

[ftalbrecht]

Copies would also make it easier to support real and imag in other implementations, e.g. those wrapping external libraries...

[sdrave]

@pmli, I have just rebased the complex branch on current pyMOR master. Could you take another look if the branch contains all the changes you need? If so, I would take one more look at performance and then merge ..

[pmli]

@sdrave I no longer have ComplexWarnings appearing, so I don't need more changes, although I would like to add a more comprehensive test.

[sdrave]

I have just merged the 'complex_merge' branch into master.