[DRAFT] Powell's algorithm

[Trombach]

Hi @stefan-k,

I'm opening this PR to seek some guidance and feedback on my first draft to implement Powell's algorithm (#222).

Currently, this code is completely untested as I am only trying to gauge whether my attempt on implementing this solver is going in the right direction.

Here are some questions I have:

1. Is my approach sensible at all? What can be improved?
2. AFAIK, Powell's method can also be implemented using Golden-section search or Brent's method. Is there a way to do this generically in argmin or will there have to be different versions of solvers for each variant?
3. I would appreciate some suggestions on a good test function to write some tests for this solver.

Thanks!

[stefan-k]

Hi @Trombach,

thanks a lot for this PR! :) At first glance it looks already quite good. Admittedly, I'm not familiar with this method. Could you provide a reference (paper, code you were following, ...) so that I can have a look, please?

1. Is my approach sensible at all? What can be improved?

I haven't seen anything that is obviously aweful, but after posting this comment I'll start a review and give some pointers where I think things could be done differently. You can also have a look at the results of the failed CI piplines, where clippy complained about a few minor things.

2. AFAIK, Powell's method can also be implemented using Golden-section search or Brent's method. Is there a way to do this generically in argmin or will there have to be different versions of solvers for each variant?

Do you mean as a replacement for the line search? It would be great to support both Golden-section search and Brent's method. One option would be to implement Linesearch for both; however I'm not sure how to implement search_direction for these. Alternatively, we could consider creating an enum with Brent(...), GoldenSectionSearch(...) and LineSearch(...) variants. Then you'd store this enum in your struct rather than the line search. Then I think we'd get away with separate constructors rather than different versions of the solver.

3. I would appreciate some suggestions on a good test function to write some tests for this solver.

I'll think about this when I'm more familiar with the method itself :)

[Trombach]

Hi @stefan-k,

thanks for the feedback! Sorry for being slow to respond, but I'm doing this as a little weekend side project :) I will go through your review suggestions as soon as I can!

I have mainly followed Wikipedia for now, but I will got through Powell's paper as well to get a more nuanced understanding. As far as I can tell it's a gradient free method where you search along provided search directions and update your list of search directions until you converge onto a minimum.

1. Thanks, I will make sure that workflows will complete in the final PR.
2. Right, an enum might make sense here. As far as I can tell any method that can find a minimum along a provided search direction without calculating gradients should work.

[stefan-k]

thanks for the feedback! Sorry for being slow to respond, but I'm doing this as a little weekend side project :) I will go through your review suggestions as soon as I can!

There's no rush, take your time :)

I have mainly followed Wikipedia for now, but I will got through Powell's paper as well to get a more nuanced understanding. As far as I can tell it's a gradient free method where you search along provided search directions and update your list of search directions until you converge onto a minimum. [...] 2. Right, an enum might make sense here. As far as I can tell any method that can find a minimum along a provided search direction without calculating gradients should work.

Thanks, I'll have a look at these. Makes sense that it works with any 1D method, but when using linesearches, it is not really a gradient free method anymore (for most line searches at least). That's fine, we just need to make that clear in the documentation :)

[stefan-k]

Hi @Trombach,

Is there any news regarding this PR? Feel free to get in touch if there is anything I can help you with :)

[Trombach]

Hi @stefan-k, My apologies for the long silence, but other things have come up. However, I'll be working on improving this PR over the weekend. I am going to respond to some of your comments individually. Good news: Powell published his paper with a test function so I will just be using the same for code tests.

[stefan-k]

My apologies for the long silence, but other things have come up.

No worries :)

Good news: Powell published his paper with a test function so I will just be using the same for code tests.

This is indeed very helpful :)

[Tastaturtaste]

Hi @stefan-k,

thanks for the feedback! Sorry for being slow to respond, but I'm doing this as a little weekend side project :) I will go through your review suggestions as soon as I can!

I have mainly followed Wikipedia for now, but I will got through Powell's paper as well to get a more nuanced understanding. As far as I can tell it's a gradient free method where you search along provided search directions and update your list of search directions until you converge onto a minimum.

1. Thanks, I will make sure that workflows will complete in the final PR.

2. Right, an enum might make sense here. As far as I can tell any method that can find a minimum along a provided search direction without calculating gradients should work.

How did you get access to the paper? I haven't found a publicly accessible version and also cannot access it through my research institution. I would like to take a look at it, too.

[Tastaturtaste]

I have some remarks regarding some of the statements regarding the line search used internally.

2. Right, an enum might make sense here. As far as I can tell any method that can find a minimum along a provided search direction without calculating gradients should work.

That should be correct. Powell's method is a gradient free optimization algorithm, meaning the linesearch used should also not need gradient information. A typical choice according to wikipedia for such line search algorithms would be Golden-section-search or Brent's method. Another important characteristic the line search method has to fullfill is that it needs to be a bi-directional line search, meaning the coefficient of the search direction can also be negative. Considering these points maybe it would make sense to hard code a specific line search method or make a new trait BidirectionalLinesearch which has to be fullfilled by the use supplied line search method. Scipy has gone the way of hard coding the Brent method for the unbounded case and (as far as I can tell) the Golden-section-search for the bounded case. Currently I am working on a direct translation of scipy's implementation based on @Trombach's work.

[stefan-k]

How did you get access to the paper? I haven't found a publicly accessible version and also cannot access it through my research institution. I would like to take a look at it, too.

Did you manage to get hold of it? I was able to access it via my institution.

That should be correct. Powell's method is a gradient free optimization algorithm, meaning the linesearch used should also not need gradient information. A typical choice according to wikipedia for such line search algorithms would be Golden-section-search or Brent's method. Another important characteristic the line search method has to fullfill is that it needs to be a bi-directional line search, meaning the coefficient of the search direction can also be negative. Considering these points maybe it would make sense to hard code a specific line search method or make a new trait BidirectionalLinesearch which has to be fullfilled by the use supplied line search method. Scipy has gone the way of hard coding the Brent method for the unbounded case and (as far as I can tell) the Golden-section-search for the bounded case.

Thanks for the clarification. I like the idea of a BidirectionalLinesearch trait. We already have Brent's method and Golden section search, therefore they could just implement the new trait.

Currently I am working on a direct translation of scipy's implementation based on @Trombach's work.

That's great! I don't know what @Trombach's plans are regarding this PR, but it seems as if this PR has gone stale. I'm therefore happy to consider another PR by you.

[Tastaturtaste]

Did you manage to get hold of it? I was able to access it via my institution.

Sadly not, but I think following the scipy implementation shouldn't be to bad for now. The algorithm can still be modified later after all.

That's great! I don't know what @Trombach's plans are regarding this PR, but it seems as if this PR has gone stale. I'm therefore happy to consider another PR by you.

I have to see how much time I can invest at the moment, but will try my best. Unfortunately I can't compile the test suite at the moment. I get the following error in a dependency named netlib-src:

   Compiling lax v0.16.0
error: could not find native static library `blas-netlib`, perhaps an -L flag is missing?

error: could not compile `netlib-src` due to previous error
warning: build failed, waiting for other jobs to finish...

I don't know yet what to do to solve this issue.

[Tastaturtaste]

Unfortunately I can't compile the test suite at the moment. I get the following error in a dependency named netlib-src:

   Compiling lax v0.16.0
error: could not find native static library `blas-netlib`, perhaps an -L flag is missing?

error: could not compile `netlib-src` due to previous error
warning: build failed, waiting for other jobs to finish...

I don't know yet what to do to solve this issue.

I found a workaround for now. It seems the issue is with a ndarray-linalg feature on windows activated with this line in the Cargo.toml from argmin-math:

_dev_linalg_0_16 = ["ndarray-linalg_0_16/netlib-static"]

The README from the ndarray-linalg repo say on windows only the Intel MKL backend is supported. Changing the cited line to

_dev_linalg_0_16 = ["ndarray-linalg_0_16/intel-mkl-static"]

as well as deleting the netlib feature from ndarray-linalg in the Cargo.toml from argmin works for me. Since this backend is supported on all platforms, maybe it should be made to be the default for everyone?

[stefan-k]

Sadly not, but I think following the scipy implementation shouldn't be to bad for now. The algorithm can still be modified later after all.

Feel free to get in touch with me via stefan.kroboth AT gmail.com. I may have a solution.

It seems the issue is with a ndarray-linalg feature on windows ...

Oh, yes, that won't work on windows. In the past, netlib was the only thing I could get to work reliably in the CI, but things may have changed for the better in the meantime. If you choose to open a PR we will see if MKL works. Anyway apologies for the inconvenience and thanks for the fix :)