# how can I do a double summation with dependent indices?

4 views (last 30 days)
Emanuele Biondini on 11 Jun 2020
Commented: Emanuele Biondini on 15 Jun 2020
Hi everyone, I have a complicated problem: I have the m and t matrix of 700 lines and I have to carry out this summation that I have to minimize. I tried to do the summation before turning the result into a function handle. However I can't do the minimization because I get the message "out of memory". Does anyone have any ideas on how to do this considering that I have to minimize the result of the summation? k, alpha, mu, c and p are the parameters that I have to find using the minimization algorithm. (N = 700). Thank you all!

Jeff Miller on 12 Jun 2020
What is that m_c in the numerator of the inner summation? Are you trying to minimize this with fminsearch or what? Perhaps it would help if you posted some of your code. I don't see anything here that would cause an out of memory error.
Emanuele Biondini on 12 Jun 2020
Attached is the matlab code that I have developed up to now. As I said I have to do the optimization of the summation obtaining the minimum point of Summ1. The values of a, c, d, k and u are the values of the unknown parameters that I am looking for. As you can see, I first carried out the double summation, then I transformed it into a function handle. But I hope there is a faster way to do it because the time for the calculation is really very long.
The double summation is:
The second problem is that when I go to do the optimization with 'simulannealbnd' matlab goes to occupy all the RAM (16Gb) while 29% of the CPU. When 99% of the RAM is occupied the PC crashes.
I therefore hope that we can do these operations in another way.
Thank you!
syms a c k d u
Ts=Ctemp(1,15);
Te=Ctemp(length(Ctemp),15); %t=Ctemp(:,15) & m=Ctemp(:,10)
Summ=0;
Su1=0;
mc=2.5;
for i=2:length(Ctemp) %length(Ctemp)=784
i
Su1=0;
for j=1:i-1
Su1=Su1+(k.*exp(a.*(Ctemp(j,10)-mc)))./(Ctemp(i,15)-Ctemp(j,15)+c).^d;
end
Summ=(Summ+log(u+Su1));
end
Summ1=-Summ;
f=matlabFunction(Summ1);
f = string(char(f));
f = replace(f,"@(a,c,d,k,u)","@(x)");
f = replace(f,["a","c","d","k","u"],["x(1)","x(2)","x(3)","x(4)","x(5)"]);
f11 = str2func(f);
%%OPTIMIZATION
x0=[0.001, 0.001,0.01,0.01,0.01]
lb=[0.001,0.001,0.01,0.01,0.01]
ub=[1000,10,1000,2,1000]
options = optimoptions('simulannealbnd','PlotFcns',...
{@saplotbestx,@saplotbestf,@saplotx,@saplotf});
x=simulannealbnd(f11,x0,lb,ub,options)
Emanuele Biondini on 12 Jun 2020
Thank you Jeff! You have been very kind to help me, your code is perfect!
In general I prefer to use "simulannealbnd" just because I can insert boundaries (upper boundary and lower boundary). For this analysis the boundaries are fundamental as the parameters cannot be <= 0 and, for example p must be between 0.3 and 2. If I carry out the minimization with "fminsearch", it gives me some values of the negative parameters and I wish it didn't happen (all parameters should be positive and within the range provided).
how can we consider the bounds?
Thank you!
ests =
4.0279 18.7059 13.9959 -0.0001 -25.1209
% x(4) and x(5) are negative but it must be positive

Jeff Miller on 12 Jun 2020
Sorry, I am having trouble relating your code to your equation, and I would have no idea how to do this with syms. But maybe good old fminsearch will help, something like this:
m = rand(10,1)*3; % I generate some random values to use for m & t.
t = rand(10,1)*3; % I guess you have real values to use.
startXs = [1 1 1 1 1]; % starting values for u,k,alpha,c,P, in that order (wild guesses)
% here is the function you want to minimize.
% fminsearch will adjust the x values
% m & t are your "data" set above
fun2Minimize = @(x) dblSumFun(x(1),x(2),x(3),x(4),x(5),m,t);
% call fminsearch and it will try to find the best values
% for u,k,alpha,c,P
ests = fminsearch(fun2Minimize,startXs)
function dblSum = dblSumFun(u,k,alpha,c,P,m,t)
% the function that you want to minimize
dblSum = 0;
mc = 2.5;
for i=1:length(m)
jSum = 0;
for j=1:i-1
jTerm = k * exp( alpha*(m(j) - mc) ) / (t(i) - t(j) + c)^P;
jSum = jSum + jTerm;
end
iTerm = log(u + jSum);
dblSum = dblSum + iTerm;
end
end
Note that I have ignored the bounds that you specified. Maybe you won't need those if the best parameters are inside the boundaries. If you do need them, a little more code will be required inside dblSumFun. But I don't want to elaborate because maybe this is not at all what you want.
By the way, don't start the parameter values near the boundary as you seem to have done in your sample code. Start at the most plausible values you can think of, which are presumably more in the middle of the range.

Emanuele Biondini on 12 Jun 2020
Thank you Jeff! You have been very kind to help me, your code is perfect!
In general I prefer to use "simulannealbnd" just because I can insert boundaries (upper boundary and lower boundary). For this analysis the boundaries are fundamental as the parameters cannot be <= 0 and, for example p must be between 0.3 and 2. If I carry out the minimization with "fminsearch", it gives me some values of the negative parameters and I wish it didn't happen (all parameters should be positive and within the range provided).
how can we consider the bounds?
Thank you!
ests =
4.0279 18.7059 13.9959 -0.0001 -25.1209
% x(4) and x(5) are negative but it must be positive
Jeff Miller on 13 Jun 2020
Well, maybe you can just use the dblSumFun with simannealbnd. But if not...
A standard trick for enforcing bounds with search routines is to use a function that maps any real number into a desired range. For example,
function trans = Real2Bounded(minimum,maximum,real)
% Convert an arbitrary real number to a number within certain bounds.
trans = (maximum - minimum) ./ (exp(-real*0.005) + 1) + minimum;
if isnan(trans)
fprintf('Real2Bounded error at min = %f, max = %f, real = %f\n',minimum,maximum,real);
end
end
One (pretty ugly) way to use this function is to change
function dblSum = dblSumFun(u,k,alpha,c,P,m,t)
to
function dblSum = dblSumFun(uReal,kReal,alphaReal,cReal,Preal,m,t)
u = Real2Bounded(0.001,1000,uReal);
k = Real2Bounded(0.001,1000,kReal); % and so on
If you do that, then you have to realize that fminsearch is not modifying your actual parameter values but instead modifying these reals that are transformed into your parameter values. To accommodat that, you have to use the reverse transform in dealing with fminsearch:
function real = Bounded2Real(minimum,maximum,trans)
% Inverse of Real2Bounded: Convert a real number within certain bounds to an arbitrary real number.
real = log( (trans-minimum) ./ (maximum-trans)) / 0.005;
if isnan(real)
fprintf('Bounded2Real error at min = %f, max = %f, trans = %f\n',minimum,maximum,trans);
end
end
For example, the starting values that you pass to fminsearch need to be transformed into the real scale that it will be using:
uRealStart = Bounded2Real(0.001,1000,uStart); % the arbitrary real for that u
startXs = [uRealStart, etc]
Then, after fminsearch runs, you need to transform its final x's back into your desired parameter estimates:
ests = fminsearch(fun2Minimize,startXs)
uFinalest = Real2Bounded(0.001,1000,ests(1)); % and so on
It's tedious to do this for all the parameters, but it works.
Emanuele Biondini on 15 Jun 2020