Fast Element-Wise power to non-integer value

I am working on a non-linear optimization problem using fmincon. My current setup delivers the best results so far in simulations.

Now I am working on optimizing my code and I have a bottle-neck in computation time in the expression:

(1-(1-x.^a).^b);

where 'x' is a vector of Nx1 observations, with N "medium-large" (> 30000 rows). Moreover, 'a' and 'b' are non-integer parameters > 0 (to be fitted)

f is the function that I am minimizing, thus, it's called by the optimizer several times (90000 in my current setting) resulting in a significant bottleneck in the overall procedure.

Are there any ways to improve the speed of this calculation? According to the code, profiler, this line is responsible for 25% of the total computation time in seconds. A replicable example is shown below:

%Example
N = 30000;
x = linspace(0,1,N)'; %vector of N x 1
a = 1.2;        %positive non-integer power
b = 0.3;        %positive non-integer power

f3 = @() (1-(1-x.^a).^b);

timeit(f3)

ans = 0.0029

Although it seems fast, if the function is called by the optimization procedure 90000 times, then the total time is around 0.0029*90000 = 261 seconds

Question:

Is there any way to make this computation faster? So far, I've tried

1. bsxfun, however it's included in the .^ command for newer releases
2. Increase maxNumCompThreads to 8
3. Relax optimization Tolerance levels
4. Tried other optimizers (such as lsqnonlin) but deliver worst estimates
5. write the expression as logarithm and then take exponential (it seems slower in preliminar trials)

Thank you for your comments

All best!

[Bruno Luong]

If you want to find 2 parameters to fit 30000 data, you could probably some sort of agregate your data and fit with a less overdertermined system, and quicker to evaluate.

You could use such strategy i a hirarchy of problems: to approach the solution, then using the solution to starting point of another lesser agregated modified problem, and make the convergence with less number of iterations than 90000.

[John D'Errico]

Big problems take big time. Why would they have implemented a slower method than possible of raising a number to a power?

Note that for large arrays, MATLAB will already automatically use multiple cores. However, 30000 is not really that large. In a quick test, I tried a test with a 20000x20000 matrix, then raised those elements to a non-integer power. It took that large of a matrix before I started to see multiple cores being seriously used.

[Alfonso Silva]

@John D'Errico @Bruno Luong

Thank you for your replies and for your time. Indeed big problems require some time. I've managed to reduce the time by avoiding unncessary calculations when possible, but ultimately it's about waiting.

Thank you again!

[Walter Roberson]

[I was wondering whether it would be effective to taylor f3 to a reasonable number of terms. I experimented with taylor around b = 0.5 (midpoint of the range), tested at 0.1 0.2 0.3 0.4... and the results for order 10 or order 20 were not encouraging.

Answer 1

[Torsten]

I don't know your optimization problem, but this helps saving function evaluations if objective function and constraints share the same function:

https://de.mathworks.com/help/optim/ug/objective-and-nonlinear-constraints-in-the-same-function.html

[John D'Errico]

Sadly, there are probably no constraints beyond bounds, if lsqnonlin was also an option.

[Torsten]

Ah, I missed that. Thank you @John D'Errico

[Alfonso Silva]

No complex constraints are included. Only bounds for parameters to ensure economic relevance

Thank you1

[Alfonso Silva]

I'll accept this answer since, the idea of avoiding unnecessary evaluations did work

Thanks!

[Walter Roberson]

Hmmm I am surprised that guide document does not refer to memoize()