feat: Add Dykstra's algorithms and tests.

[tttamaki]

I added the following features. The large part is in three files, and a separate test file.

• DykstrasProjection in projection/DykstrasProjection.py
  • (Parallel) Dykstra's projection algorithm computing convex projection $P_C(x)$ of $x$ to the intersection of convex sets $C = \cap_{i=1}^m C_i$
• DykstrasProjectionProx in proximal/DykstraLikeProximal.py
  • The corresponding indicator function (prox)
• DykstraLikeProximal in proximal/DykstrasProjectionProx.py
  • (Parallel) Dykstra-like proximal algorithm computing prox of the sum of convex functions $prox_{\tau \ f + g}$ or $prox_{\tau \ \sum_{i=1}^m w_i f_i}$
• pytests/test_dykstra.py
  • tests in a separate test file.

These are a bit large files, so I would appreciate it if you have time to review the features.

Notes

• Naming Conventions
  • I think that the class names may not follow the pyproximal convention. Please suggest better names if you have recommendations.
• Dykstra's projection
  • I have created a gist demonstrating how it works. The implementation appears to work correctly.
• Dykstra-like proximal
  • I' aware of issue Combining proximal operators #116 for Dykstra-like proximal operators, but it is still open for more than two years and does't include Dykstra's projection.
  • As noted in my code comments, I added self.use_original_tau as a hidden feature to align the code with the known equation in the literature, which doesn't pass the tests. I need more time to understand the proofs in the literature. Anyway, the modified implementation passes all tests and performs as expected. (I have only found one other library pylearn-parsimony which uses the known equation, but it lacks tests)
• About test
  • As test functions are too many and hard to be related to other tests, so I separated tests for this from others.
  • Seeds are fixed to 10 like as other tests, but I have verified that tests pass with seeds ranging from 0 to 1000 during debugging and convergence testing.

[mrava87]

@tttamaki this is awesome!

Give me a week or so to review them 😄 in the meantime, give that this is a completely new feature, would you mind adding an example or tutorial (where you see more fit - our rule of thumb is, if the script just applies the operator to some dummy data - a bit like a test - it’s more an example, if the script takes an interesting problem and solves it, then it’s a tutorial) so users will be able to see how to use the operators in action

[mrava87]

Ok I spoke too fast without reading the entire message… your gist would be a perfect example 😀

[mrava87]

@tttamaki I started reviewing your PR and I think these Dykstra algorithms look like a great addition to pyproximal 😄

Regarding the Github issue you mention, I wouldn't worry about too much as, like you say, it is 2 years old so I doubt anything will ever happen there (once we merge this PR I may close it).

To the code: I think in general it is already very well written, I let some comments here and there (mostly typos and unclear sentences in docstrings). However, the more I look at it, the more I think some changes may be needed to make sure that this feels more like a natural addition to PyProximal that follows the various standards that we established with our operators and solvers. In general, I feel you have kind of mixed up the purspose of these operators (which to me is Intersection of two or more projections or sum or two or more proximable function) with the algorithm used to find a solution of the underlying optimization problems. I would prefer we call those operators something like Intersection for the projections/indicator functions and maybe Sum for the generalization to any proximal function. This way the operator would give a clear idea to what they are meant to do in terms of defining an optimization problem (and its objective function) and the fact that the Dykstra algorithm is used to solve such problems is to me an algorithmic choice that we adopt but should not be exposed so explicitly to the user.

So in other terms I would like to see someone do:

Op = 
l2 = L2(Op, d)
circle = EuclideanBallProj(np.array([0.0, 3.5]), 5)
box = BoxProj(np.array([-5.0, -2.5]), np.array([5.0, 2.5]))
constr = Intersection([circle, box])

ximv = ADMM(l2, constr, ...)

more than constr = DykstrasProjectionProx([circle, box])

Another thing that I am not so convinced is why we need DykstraLikeProximal and DykstrasProjectionProx. Unless I do not understand something, _dykstra_projection in DykstrasProjectionProx is nearly identical to _dykstra_like_proximal_algorithm in DykstraLikeProximal with the only change that instead of calling self.projections[0] we call self.ops[0].prox.... could it not be possible to create a single solver (this could be in a private file called _dykstra or something) and have a flag that says if we pass proximals or projections (or even discover it from the class type as all prox operators are subclassed of ProxOperator), based on which we can do something likes this (this is just a pseudocode to give an idea):

if prox:
  self.proj_or_prox = projections
else:
  self.proj_or_prox = [partial(projection.prox, tau) for projection in projections]

Finally, one last point that we could test later, once the code is in place: what about if we add __add__ to ProxOperator to allow one doing something like this:

circle = EuclideanBall(np.array([0.0, 3.5]), 5)
circle1 = EuclideanBall(np.array([2.0, 3.5]), 3)
box = Box(np.array([-5.0, -2.5]), np.array([5.0, 2.5]))
constr = circle + circle1 + box

and internally the __add__ method will just form and return a DykstraLikeProximal/Sum operator.

What do you think?

[tttamaki]

Thank you for your review, I’m afraid that it might have taken a lot of time.

• About the naming of Intersection and Sum
  I think it is good. However, when I first considered this PR, I noticed that pyproximal.Intersection already existed, which is the name I wanted to use. (By the way, the naming of pyproximal.Intersection is confusing. When I first found it, I misunderstood it would run Dykstra’s Projection for m >= 2.)

• Why both DykstraLikeProximal and DykstrasProjectionProx exist
  You’re right. Both are the same for projections, so DykstrasProjectionProx could be unnecessary. The reasons I included both are:

  • I created both to keep that consistency as every projection function has a corresponding prox. Now, It think that having both might even help users. When someone wants a “projection to intersection,” it may feel unintuitive to call a function named “proximal.”
  • For m=2 the algorithms are the same, but for m >= 2 there is a subtle difference. In parallel Dykstra's projection, the m projections are applied sequentially. In parallel Dykstra-like prox, the m prox results are averaged and redistributed. Both are proven to converge, so practically the results should be the same. I think that having both might be useful for debugging, e.g. checking if the projections match.

  If these reasons aren't compelling, I agree to remove DykstrasProjectionProx.

• About Sum and __add__
  I agree that Sum is clear as a name. It should take a list of proximable functions as arguments.
  On the other hand, the idea of __add__ is probably not good because it is inefficient. My understanding is that operators like __add__ can only be implemented with two operands at a time. So while this syntax is easy to read:

  constr = circle + circle1 + box

  internally it becomes:

  constr = Sum(Sum(circle, circle1), box)

  with the version with m=2. This means the outer Sum calls the inner Sum at every iteration, leading to nested loops. As the number of terms grows, the cost can grow exponentially.

  In contrast, if Sum accepts a list:

  constr = Sum([circle, circle1, box])

  then the parallel version (m >= 2) runs and only a single loop is required, so the cost grows linearly with the number of terms.

[mrava87]

@tttamaki oh yeah I forgot about Intersection already existed. Well it is an intersection but a special one as it is only for a given type of convex set. It comes form Appendix A in https://cvg.cit.tum.de/_media/spezial/bib/chambolle_et_al_siims12.pdf and I had implemented when I was working on a segmentation project as it was suggested to be a good constraint for segmentation. It runs the Djkstra algorithm under the hood but I think it’s a bit of a special version of it tailored to the specific problem - the paper explains this in more details if interested.

So, based on the above, and the fact we can’t really change the name of the current operator called Intersection, what about calling your one GenericIntersection, where Generic is used to indicate that it can work with any convex set whose projection is implemented in one of our Projection operator?

To the second point about why both ‘need DykstraLikeProximal and DykstrasProjectionProx exist’ (note that what I asked doesn’t match with what you replied)… yes I am keen to keep consistency so to still have a projection operator and then a proximal for its indicator function that just calls the projection but what I’m not sure is about DUPLICATING the actual implementation of the Djkstra algorithm. So in other words I am fine to have GenericIntersection and Sum in the proximal submodule but internally they should both call the same algorithms (as I said about, it seem to me their algorithms are identical other than one is calling the call method of a projection and the other the prox method of a proximal). The actual Djkstra algorithm could but put into a _Djkstra file in optimization if we don’t want to expose it directly or in primal if we think it may be also worth to allow one to just run the algorithm passing some projection/proximal ops, and then internally we call this solver…

Ok, good to hear that I understood correctly. For m>=2 the algorithms are slightly different. Fine, if this is what is in the literature we can definitely do that, but I don’t see an issue as we can always implement all variations in _djkstra and the call the appropriate one in the different operators.

For the add, I see your point and I had initially the same thought. But, although I will need to play a bit with it to see if I’m right, my idea was that in the add we can easily check if one of the operators is already of type sum, extract its terms from self.ops and re-initialize a new Sum with all previous ops plus the new one… so ultimately what is returned is a single Sum and the prox operator will run a single Djkstra algorithm 😀

Hope this makes sense and you can start making the suggested changes? Ping me when you want me to do a second review 😉

[tttamaki]

@mrava87 Thank you for the clarification. I'll try to modify the PR. But the __add__ looks a bit tricky, so give me time for while.

[mrava87]

@tttamaki no worries, there is no rush 😀 for the add, feel free to have a go at it, but if you prefer I can do it once the rest of the PR is ready

[tttamaki]

@mrava87 Thank you. I have finished to refine the PR except the add, so please check if I understand your suggestions correctly.
For the add, I found that ProximalOperator has already __add__ for a particular purpose, so I think that seamlessly integrating the idea to the current implemeantation would be horribly tricky. So I think modifying the add is better to be a separate task.

[mrava87]

@tttamaki this looks great!

I just made a few minor changes here and there to be more consistent with the rest of the codebase (nothing in the code itself, mostly docstrings and examples).

For the __add__, let me take it in a separate PR. What I am thinking is that inside the current __add__ we can easily add logic based on the dytpe of the inputs to decide what to do - if np.array dispatch to affine addition (what is done now), if both are prox operator call the Sum passing the two prox operators, if any of them is a prox operator and the other is already a Sum, unpack the prox operators, put them into a single list and pass them to Sum again 😄

[mrava87]

@tttamaki I had a look at wrapping Sum into __add__ but whilst it seems doable I am not so convinced anymore whether it is useful, as Sum takes quite a few extra parameters and those would need to be kept default (or taken from the current Sum if you do ProxOperator+Sum / Sum+ProxOperator)... I'll give it another thought and maybe try a bit more to see if it makes sense or not to have such convenience..