Projecting complex LincombOperator onto range_basis

[alexandre-pasco]

Describe the bug When projecting a LincombOperator with complex coefficients onto some range_basis using pymor.algorithms.projection.project, the result is not the one expected.

To Reproduce

from pymor.operators.numpy import NumpyMatrixOperator
from pymor.operators.constructions import LincombOperator 
from pymor.algorithms.projection import project
from pymor.parameters.functionals import ProjectionParameterFunctional, Mu
import numpy as np

mat =np.array([[1 + 1.j], [1]])
operators = [NumpyMatrixOperator(mat), NumpyMatrixOperator(mat)]
coefficients = [ProjectionParameterFunctional('mu', size=2, index=i) for i in range(2)]
op = LincombOperator(operators, coefficients)
mu = Mu({'mu':[1, 1.j]})
u = op.range.from_numpy(np.array([[1, 1]]))

m1 = project(op, u, None).assemble(mu).matrix
m2 = project(op.assemble(mu), u, None).matrix

print(f"Projecting then assembleing = {m1}")
print(f"Assembling then projecting = {m2}")

Expected behavior I would expect project and `assemble' to commute (in some sense).

System information:

pyMOR Version 2022.2.1

Python 3.11.3 on Linux-5.15.0-73-generic-x86_64-with-glibc2.31

External Packages
--------------------
DEALII:      missing
DUNEGDT:     missing
FENICS:      missing
GL:          missing
IPYTHON:     missing
IPYWIDGETS:  missing
MATPLOTLIB:  3.7.1
MESHIO:      missing
MPI:         missing
NGSOLVE:     missing
NUMPY:       1.24.3
PYMESS:      missing
PYTEST:      missing
PYTHREEJS:   missing
QT:          missing
QTOPENGL:    missing
SCIKIT_FEM:  missing
SCIPY:       1.10.1
SCIPY_LSMR:  present
SLYCOT:      missing
SPHINX:      missing
TORCH:       missing
TYPER:       0.9.0
VTKIO:       missing

Additional context in pymor.algorithms.projection, I think a .conj() is missing. In particular, I think the line 122 should be

return NumpyMatrixOperator(V.to_numpy().conj(), source_id=op.source.id, name=op.name)

instead of

return NumpyMatrixOperator(V.to_numpy(), source_id=op.source.id, name=op.name)

Indeed, the problem seems to be fixed when I insert a rule with the corresponding correction into ProjectRules.

[sdrave]

Thanks for the bug report @alexandre-pasco. I believe the issue is in

https://github.com/pymor/pymor/blob/84f7d613aab4380aacae08a27b9a3da461f55f13/src/pymor/algorithms/projection.py#L127

which should be

V = op.apply_adjoint(range_basis).conj()

I'll fix this as soon as I have checked that I am not overlooking anything.

[alexandre-pasco]

I think that this would fix the problem, however I also think that it would lead to another problem on line 135, return VectorArrayOperator(V, adjoint=True, name=op.name), which already take the adjoint operator. Thus, I think the fix should be on line 132, which would become return NumpyMatrixOperator(V.to_numpy().conj(), source_id=op.source.id, name=op.name) instead of return NumpyMatrixOperator(V.to_numpy(), source_id=op.source.id, name=op.name).