DFT-based Operators

[artpelling]

This PR adds a new module operators.dft for matrix-free operators with DFT based matrix-vector multiplications, namely CirculantOperator, HankelOperator and ToeplitzOperator. See here to get the gist.

In contrast to NumpyHankelOperator the operators can now also be non-square.

Todo:

• more tests
• finish the implementation of the block versions (to_matrix rules and tests)
• make ERA work with the new module

[artpelling]

I took some time today to polish my code a bit. We've talked about this a while ago. I decided to call them dft-based operators in order to avoid confusion with the "structured operators" that may come.
Let me know what you think about it!
I am planning to add more system identification methods before/at pyMOR School that will make more use of these operators.

[pmli]

It seems ok to me so far. It would be good to resolve conflicts.

[artpelling]

It seems ok to me so far. It would be good to resolve conflicts.

Cool. Yeah, I forgot to finish the rebase.

[artpelling]

I've decided to remove BlockCirculantOperator, BlockHankelOperator and BlockToeplitzOperator because they did not really add anything important other than a convenience function for construction. I think in these rare cases, it is not hard to construct them by hand.

I've finalized BlockDFTBasedOperator with tests and a to_matrix rule.

The functionality right now is exactly what I need, but in regards of code structure I am not happy yet.
The operator is very similar to BlockOperator, however, it does not get blocks and glues them together to a block matrix. It is constructed by a 2d-array of Operators that are then interlaced, e.g. $A,B,C,D\in\mathbb{C}^{2\times2}$ get turned into

$$\begin{bmatrix} a_{1,1}&b_{1,1}&a_{1,2}&b_{1,2} \\\ c_{1,1}&d_{1,1}&c_{1,2}&d_{1,2} \\\ a_{2,1}&b_{2,1}&a_{2,2}&b_{2,2} \\\ c_{2,1}&d_{2,1}&c_{2,2}&d_{2,2} \end{bmatrix}$$

instead of

$$\begin{bmatrix} A&B\\C&D \end{bmatrix}.$$

The nice thing is, that $A,B,C,D$ can be arbitrary operators and each operator's apply method is called in BlockDFTBasedOperator making it efficient if some of theses operators are matrix-free. Furthermore, the source and range dimensions of the individual operators do not have to match. In this case, to_matrix and apply fill missing entries with zeros and minimal computation is done.

Right now, BlockDFTBasedOperator inherits from DFTBasedOperator which does not make sense. How could I call this Operator and where does it fit? I would think operators.block but that could cause confusion. Any ideas @pymor/main ?

[pmli]

I'm thinking that maybe these Operators should be named Numpy... for consistency since the source and range are NumpyVectorSpaces.

In that sense, BlockDFTBasedOperator would not fit in pymor.operators.block if the source and range are not BlockVectorSpaces.

[artpelling]

I'm thinking that maybe these Operators should be named Numpy... for consistency since the source and range are NumpyVectorSpaces.

Thats true. I'm happy to change the names, maybe even of the module.

In that sense, BlockDFTBasedOperator would not fit in pymor.operators.block if the source and range are not BlockVectorSpaces.

Yes. But what can we call it? It still does a lot of things similarly. It might be possible to route the operations in apply through a BlockOperator by apropriately mangling and summing input and output vectors... But I don't really like that.

[pmli]

In that sense, BlockDFTBasedOperator would not fit in pymor.operators.block if the source and range are not BlockVectorSpaces.

Yes. But what can we call it? It still does a lot of things similarly. It might be possible to route the operations in apply through a BlockOperator by apropriately mangling and summing input and output vectors... But I don't really like that.

Would it make sense to have a single Operator for both scalar and blocked versions? E.g., HankelOperator that accepts both a scalar or a matrix sequence?

[artpelling]

Would it make sense to have a single Operator for both scalar and blocked versions? E.g., HankelOperator that accepts both a scalar or a matrix sequence?

That's how it was before. The old NumpyHankelOperator accepted one or three-dimensional arrays. This new version is more flexible: One can mix arbitrary Operators, so especially HankelOperators with different "lengths". In apply only DFTs of necesssary length are performed for each "channel". With the old approach, one would have to zero pad all "channels" to a common length. Similar arguments can be made for to_matrix.
Furthermore, it might be handy in the future to have some "channels" as ZeroOperators etc.

A compromise could be to bring back NumpyBlockHankelOperator with the from_numpy classmethod and keep NumpyBlockDFTBasedOperator as parent thats not intended for the user..

[pmli]

It looks like the new Operators are gone.

[pmli]

Concerning the non-blocked Operators, it seems to me that composition would be better than inheritance. Since both Toeplitz and Hankel matrices use circulant matvec, my suggestion would be to merge NumpyDFTBasedOperator and NumpyCirculantOperator into a single NumpyCirculatOperator class that inherits from Operator and CachableObject. Then NumpyToeplitzOperator and NumpyHankelOperator would inherit from Operator and instantiate a NumpyCirculantOperator in their __init__ methods which would be used in apply.

I'll think about the blocked Operators more later.

[artpelling]

Concerning the non-blocked Operators, it seems to me that composition would be better than inheritance. Since both Toeplitz and Hankel matrices use circulant matvec, my suggestion would be to merge NumpyDFTBasedOperator and NumpyCirculantOperator into a single NumpyCirculatOperator class that inherits from Operator and CachableObject. Then NumpyToeplitzOperator and NumpyHankelOperator would inherit from Operator and instantiate a NumpyCirculantOperator in their __init__ methods which would be used in apply.

You are right. Mathematically speaking, a circulant matrix is a special kind of Toeplitz matrix, so I tried to reflect that in the inheritance structure. Probably not important though...

I'll think about the blocked Operators more later.

I am currently trying to run the code on a huge example and it turns out that my approach is very inefficient. I think the speed benefits of using 2d-FFTs on large arrays is much better than avoiding unecessary zeroes. So I am working towards your suggestion (close to the old code) where e.g. HankelOperator accepts one and three-dimensional arrays.

[artpelling]

@pmli I've refactored the code and added tests also to fixtures/operators.py. I think this version is fast and does not use a lot of code. I've improved the matvec loop in NumpyHankelOperator and used scipy.fft instead of numpy.fft. This enables the use of the scipy.fft.set_workers context manager which can be used for parallelization.

[sdrave]

@artpelling, there seems to be an issue with the api docs (see ci build). Can you take a look, or do you need help?

[artpelling]

@artpelling, there seems to be an issue with the api docs (see ci build). Can you take a look, or do you need help?

Fixed it. Also had to fix some operator fixtures. They revealed a minor but fatal mistake in the apply method that I have fixed as well :)

[artpelling]

@sdrave Do you have an idea what is wrong with the tests? almost_equal fails for the new operators in test_InverseOperator and test_InverseAdjointOperator although both results seem to be identical. When I mimic the tests in a separate script, they are identical. When I add a print statement in the tests, the tests pass...

[pmli]

I think I get the issue. For an Operator that doesn't implement apply_inverse and to_matrix works, the first call to apply_inverse will use NumpyMatrixOperator(to_matrix(...)).apply_inverse, and all later calls will use lgmres. This makes the tests fail when the two solutions are different enough, but solving again will produce exactly the same solution.

In Operator.apply_inverse:

pymor/src/pymor/operators/interface.py

Lines 236 to 263 in efd4095

	if self.linear:
	if solver_type is None or solver_type == 'to_matrix':
	mat_op = None
	if not hasattr(self, '_mat_op'):
	if solver_type is None:
	self.logger.warning(f'No specialized linear solver available for {self}.')
	self.logger.warning('Trying to solve by converting to NumPy/SciPy matrix.')
	from pymor.algorithms.rules import NoMatchingRuleError
	try:
	from pymor.algorithms.to_matrix import to_matrix
	from pymor.operators.numpy import NumpyMatrixOperator
	mat = to_matrix(assembled_op, mu=mu)
	mat_op = NumpyMatrixOperator(mat)
	if not self.parametric:
	self._mat_op = mat_op
	except (NoMatchingRuleError, NotImplementedError):
	if solver_type == 'to_matrix':
	raise InversionError
	else:
	self.logger.warning('Failed.')
	if mat_op is not None:
	v = mat_op.range.from_numpy(V.to_numpy())
	i = None if initial_guess is None else mat_op.source.from_numpy(initial_guess.to_numpy())
	u = mat_op.apply_inverse(v, initial_guess=i, least_squares=least_squares)
	return self.source.from_numpy(u.to_numpy())
	self.logger.warning('Solving with unpreconditioned iterative solver.')
	return genericsolvers.apply_inverse(assembled_op, V, initial_guess=initial_guess,
	options=options, least_squares=least_squares)

the issue is that the if not hasattr(self, '_mat_op'): is missing an else branch that would set mat_op to self._mat_op.

The question then is, is it ok for apply_inverse to call to_matrix in the case of DFT Operators?

[artpelling]

This should be ready if tests would not fail. I have no idea what is happening, is this related to artifacts? @sdrave @pmli

[artpelling]

This should be ready if tests would not fail. I have no idea what is happening, is this related to artifacts? @sdrave @pmli

NVM, now nothing fails :)

[sdrave]

I haven't looked into the mathematical details. But structurally, the code looks good!