Change all linear solvers to accept single operator or list of operators

[adler-j]

Currently some solvers (e.g. pdhg) only support a single operator, while some (e.g. douglas_rachford_pd) only support list of operators.

Converting between these is trivial but takes a few lines of code that becomes repetetive. How about we add functionality to allow users to feed either of these? Converting inside ODL should be rather easy.

[aringh]

I am all for harmonizing the interfaces!

One thing that could be done is simply to explicitly extend the pdhg algorithm to the same kind of problems that is solved by douglas_rachford_pd. In that way people would not have to create product space operators for solving problems with multiple operators when using pdhg.

[adler-j]

Indeed that is half of this issue :) The other half is simply to teach douglas_rachford_pd to regard operator as [operator], and we're done! (conditions apply, e.g. we also need to fix step sizes etc)

[kohr-h]

Do we really need two interfaces with the same functionality? We saw with the vector step sizes that they add quite a bit of complexity, just to make PDHG work the same as Douglas-Rachford regarding steps.

I would advocate doing the op -> [op] thing but retiring the vector step sizes in favor of a list of operators with a list of step sizes.

[adler-j]

Frankly I'm not 100% up to date on vector step sizes, I'll have to dig deeper into that. Is the worry simply growth of complexity?

My proposal is that we make the "behind the scenes" code simply work on lists of operators in all cases, then we should only need code that converts from single op to list of op in the start.

[sbanert]

I would advocate doing the op -> [op] thing but retiring the vector step sizes in favor of a list of operators with a list of step sizes.

Now I’m triggered. :wink: Giving up the vector-valued stepsizes would make certain preconditioning tedious, intransparent and less optimised. If one would e.g. use the stepsize [1, 2] on an L¹ functional on ℝ², one would have to split ℝ² into ℝ × ℝ and use L¹ (i.e., absolute value) on both of them and a solver would have to traverse a list instead of doing the numpy pointwise multiplications.

[kohr-h]

Hehe, good. I was hoping to get some reaction. But then, if we have lists of operators and lists of step sizes, why shouldn't it be allowed that some of them are vectors?

[sbanert]

Yes, it should, but on the other hand, it is a bit of an effort to implement this freedom: Imagine a SeparableSum functional on a ProductSpace consisting of an ℝ² and an ℝ³ factor. Then we have the following choices of a stepsize given to the proximal of the SeparableSum functional:

• one float,
• a space element of the ProductSpace,
• a list of floats,
• a list of a float and an element of ℝ³,
• a list of an element of ℝ² and a float,
• a list of an element of ℝ² and an element of ℝ³ (which has the same meaning as an element of the ProductSpace).

In the implementation, this has to be somehow considered, but for this we need more than the simple np.isscalar, I guess.

Edit: I was mixing up ProductSpace and SeparableSum.

[mehrhardt]

I am a bit lost. Why would Douglas-Rochford have a list of operators as input? Isn't this usually formulated with one operator only? Of course this operator can be a ProductSpaceOperator but why do we ask for this explicitly? There are algorithms like SPDHG where it does not make any sense if the operator is not a ProductSpaceOperator but why do we now start using "lists" for these?

[adler-j]

Well I guess it's for historical reasons. The article we followed, used a list of operators. Overall though it's easier if all of the solvers support the same syntax.

For some things, e.g. TV-reg, having to create the ProductSpaceOperator wastes quite a few lines of boilerplate code. In those cases allowing users to simply give a list of operator help.

[mehrhardt]

Well, once you have a list you can easily make this a BroadcastOperator, no? Wouldn't change the number of lines too much but makes the code much simpler. Having the facility of odl, I would use it and in any case I would not make this the "default".

All of this being said, I would not mind odl operators A supporting A[0] and len(A) which should return the operator itself for a normal operator and a "component" if A is a ProductSpaceOperator. I.e. every operator is a (trivial) ProductSpaceOperator. Then the two notions discussed here are kind of the same anyway.