Speeding up the fminbnd function

Hi!

I have a function of T as,

CV_inf = @(T) T.*(4/(n+1)) - integral(@(t) psitildasq(t),(-T),(T),'Arrayvalued',true);

where

psitildasq = @(t) (1/n^2)*(sum(cos(x*t))).^2 + (1/n^2)*(sum(sin(x*t))).^2;

x is a (1*n) vector of inputs.

I want to find the global minimum as well as the first(smallest) local minimum of this function.

To find the global minimum I use

[Tinf,Tinf_val,Tinf_exitflag,output] = fminsearch(CV_inf,0);

and to find the local minimum I follow this method,

[Tloc1,Tloc1_val,Tloc1_exitflag] = fminbnd(CV_inf,0,Tinf);

    g = Tinf - Tloc1;

    if g <= 0.0001
        Tloc = Tinf;
        Tloc_val = Tloc1_val;
        Tloc_exitflag = Tloc1_exitflag;
    else
        while g > 0.0001
            [Tloc,Tloc_val,Tloc_exitflag] = fminbnd(CV_inf,0,Tloc1);
            g = Tloc1 - Tloc;
            Tloc1 = Tloc;
        end
    end

But when I have a large dataset like n=10000, the code seems to run forever. Could there be a possibility of speeding up the process?

Thanks.

Answer 1

[Walter Roberson]

If you factor the 1/n^2 out of the integral, then what is left can be integrated into a closed form that has a predictable pattern. You can calculate based upon that pattern and in so doing save all of the integrations.

The indefinite integral of

sum(cos(x[k]*t), k = 1 .. n)^2 + sum(sin(x[k]*t), k = 1 .. n)^2

is (n*t + 2 * R)

where R is

sum( sum( sin((x[p]-x[q])*t)/(x[p]-x[q]), q = 2..n), p = 1..n-1 )

so sin() of the pairwise difference in x, times t, divided by the difference in x, where the first index of the subtraction is less than the second index of the subtraction.

Another way of writing 2*R would be

sum( sum( sin((x[p] - x[q])*t) / (x[p] - x[q]), q = 1 .. n), p = 1 .. n)

Notice that the differences can be calculated ahead of time. Using the full form as being easier to write out,

D = bsxfun(@minus, x(:), x(:).');

and then for any one t, sum(sum(sin(D*t)./D)). Optimize slightly by going for vector form, calculate ahead of time,

D = reshape( bsxfun(@minus, x(:), x(:).'), [], 1);  %column vector

and for any one t, sum(sin(D*t)./D)

See below for a correction

With regards to the definite integration over -T to +T, notice that since sin(-x) = -sin(x), evaluating the sum at -T is going to be the negative of the sum evaluated at T. In just the same way, the leading term n*t at t=-T would be the negative of n*t at t=T. The indefinite integral is symmetric around therefore symmetric around 0, and so going over -T to +T is going to be 2 * going over +T by itself. The definite integral is therefore going to be

2*n*T + 2 * (the sum from above, evaluated at t=T)

With these optimizations, execution should be comparatively fast provided that you can store all of D in memory.

Remember the constant multiplier of 1/n^2 that got moved out of the integral.

Note: if it is possible that two of the variables might be the same then adjustments need to be made to the integral. The term sin((p-q)*t)/(p-q) would become 0 / 0 when p = q. This value effectively has to be treated as being t in order to get the calculations to work out properly. Indeed, this is the source of the n*t term, since there would be n places where the pairwise variables would be equal. This points out a problem in my writing D in terms of bsxfun: it is going to generate 0 terms. This leads to a slight reformulation: continue to have

D = reshape( bsxfun(@minus, x(:), x(:).'), [], 1);  %column vector

but now omit the n*t leading term, and instead use

Dr = D(D ~= 0);
num_zero_pairs = n^2 - length(Dr);

and then the indefinite integral is

num_zero_pairs * t + sum(sin(Dr *t) ./ Dr)

and the definite over -T to +T would be

2 * (num_zero_pairs * T + sum(sin(Dr.*t) ./ Dr)))

again remember the 1/n^2 multiplier

for n = 10000 then length(D) would be 50005000, 50 meg entries, 400 megabytes. Do-able. I guess it would be a fair question as to whether this exact integration would possibly be more expensive than numeric integration. It might make it easier to reason about CV_inf though.

[Torsten]

If the OP wants to minimize CV_inf, he should differentiate it with respect to T and solve CV_inf'(T)=0 (using fzero, e.g.). This way, no integration is required at all.

Best wishes

Torsten.

[Anjali]

Thanks for your answers. I will try both ways.

[Anjali]

Hi! Shouldn't we extract only the upper triangular elements in the square matrix that has the differences?

D = bsxfun(@minus, x(:), x(:).');

[Walter Roberson]

You could, in which case you would need to double the results except for the diagonal. If you just use the entire matrix the value will be counted once for each of the times and the sign will work out because the difference in sign from the trig part will get cancelled by the difference in sign of the division.

Answer 2

[Torsten]

You could take the first derivative of your function CV_inf with respect to T, set the derivative to 0 and solve for T.

That way no integration is required in your calculation.

Best wishes

Torsten.

[Walter Roberson]

You will, though, need to solve for all of the 0's and use the derivative to classify them. As they are trig functions, there is a risk of an infinite number of roots, though one might be able to establish a period for them. The differences between x values act as phase shifts, so as long as the x values can all be expressed as rational numbers, you can find a LCM (Least Common Multiple) of the numbers, multiply by 2*Pi, and that should be the period for the entire psitildasq. But then the additive T of CV_inf ... ah yes, as you can put an upper bound on the sum of n cosines (or n sines), you can put bounds on the range of psitildasq... muse, muse, muse

[Torsten]

I don't think a minimizer works more efficient on a function with several local minima than a root finder on its derivative.

Best wishes

Torsten.