BUG: minimize - COBYQA is very inefficient w.r.t time.

[andyfaff]

Describe your issue.

CI tests revealed that various COBYQA tests were taking a long time. Local testing reveals that there is a lot of overhead in COBYQA internals compared to the time taken to call the objective function. I believe that there are inefficiencies in the implementation that need to be discovered and reduced.
The example below shows that only 0.2% of time is spent in the objective function for COBYQA, whereas for nelder-mead it's 33%. COBYQA gets to the solution in a smaller number of nfev (nelder-mead doesn't actually get to a good minimum), but it takes two orders of magnitude more time.

@ragonneau

Reproducing Code Example

import numpy as np
from scipy.optimize import rosen, differential_evolution, minimize
bounds = [(0, 10)] * 10
rng = np.random.default_rng(1)
x0 = rng.random(size=10)*10

start = time.time()
res =  minimize(rosen, x0, method='cobyqa')
finish = time.time()
sftime = time.time()
_ = [rosen(res.x) for i in range(res.nfev)]
fftime = time.time()

print((fftime - sftime)/(finish - start), res.nfev, finish-start)

start = time.time()
res =  minimize(rosen, x0, method='nelder-mead')
finish = time.time()
sftime = time.time()
_ = [rosen(res.x) for i in range(res.nfev)]
fftime = time.time()

print((fftime - sftime)/(finish - start), res.nfev, finish-start)

Error message

0.0019886230044629835 1087 2.6466002464294434
0.30101082386617767 1898 0.02931809425354004

SciPy/NumPy/Python version and system information

N/A

[andyfaff]

e.g. A long time is spent in Framework.get_best_index. There is a loop in that method that calls Framework.merit repeatedly. As part of merit the constraint violations are calculated.
In a subsequent line in the loop the max constraint value is calculated. This maxcv calculation reevaluates the constraint violations again. i.e. duplication of calculation.

Furthermore, it's not clear why those constraint violations even need to be calculated for a system where there are no constraints or bounds.

[andyfaff]

Most of this loop could be vectorised.

[andyfaff]

Vectorize this line

[ilayn]

Norm accepts axis keyword. No need to loop there

[nickodell]

Hi, I'm interested in helping out on this issue. Are you interested in PRs from anyone on this topic, or just the original author?

[andyfaff]

I've started a branch at https://github.com/andyfaff/cobyqa/tree/efficiency that's looking at various efficiency things. This is based on a jupyter notebook where I was line profiling things, https://gist.github.com/andyfaff/4cfdc554f9f520fb1a6a6f983025996f.

I'm yet to find the big prize.

It would be good to have the authors on board here, as the originators of the code they're familiar with it, and it's going to be easier for them. You're definitely more than welcome to contribute here, but it might be a slow slog through profiling country.

• There are a lot of expensive scipy.linalg.eigh calls.
• It's annoying that the code wraps optimize.LinearConstraint, optimize.NonLinearConstraint with problem.LinearConstraints and problem.NonLinearConstraints instead of using the original classes directly.
• I suspect that there's a lot of linear algebra going on that doesn't need to be. e.g. if there are no constraints or bounds, then the algorithm should be designed to quickly bypass all calls to evaluation of those constraints, and the linear algebra involved with processing them to direct the minimisation.
• I suspect there's a bunch of caching that could be done, that isn't.

[nickodell]

There are a lot of expensive scipy.linalg.eigh calls.

I found it was possible to reduce the time spent in eigh() by 30% by passing driver="evd" to eigh().

E.g.

eig_values, eig_vectors = eigh(a, driver="evd", check_finite=False)

Unfortunately, it also increased nfev by 3%, which may not be worth it.

I also tried passing overwrite_a=True, but I did not find any speed-up from this.

[dschmitz89]

This looks very helpful, I had a similar impression when I had a glance over COBYQA's source code. For the sake of repo hygiene, could we move the discussion to the COBYQA repo though? It must be fixed there anyway.

[mdhaber]

Please reassess the fail_slow marks (gh-20699, gh-20741) when this is resolved.

[zaikunzhang]

Hello @andyfaff and all,

First of all, thank you very much for bringing this issue up and for proposing improvements (see also cobyqa/cobyqa#152). They are extremely important to us. We cherish your suggestions and inputs.

Meanwhile, I would like to mention that COBYQA (PRIMA as well) is primarily designed for "derivative-free optimization" (DFO), aka "black-box optimization" or "simulation-based optimization". We assume that function evaluations predominate the cost. In practice, each function evaluation can take seconds, minutes, hours, or days, because it often involves sophisticated simulations. In contrast, the costs incurred inside the solver (e.g., floating point operations, memory allocations) are assumed to be negligible. That is why we commonly use the number of function evaluations to measure the performance of DFO solvers.

However, we (Tom @ragonneau and I) always keep in mind that our assumptions can be wrong, and users may use our DFO solvers to solve non-DFO problems. Therefore, we should not ignore the internal costs of our solvers. Prioritization is key:

Reducing function evaluations >> enhancing code simplicity & clarity >> economizing solver's internal costs.

Any improvement of the code that does not increase the number of function evaluations should be integrated, just like what we have done at cobyqa/cobyqa#152 (credit goes to @andyfaff . Thank you!).

If you are interested in this topic, the following notes on DFO may be interesting to you. You are invited to make comments there:

https://github.com/orgs/libprima/discussions/145

Many thanks again!

Cheers,
Zaikun