is it possible to use a lookup table for optimization?

I am doing an optimization and need to include a look up table as constraint. A variable is another's lookup table value. Is there a way I can do that?

Answer 1

[Walter Roberson]

Beyond what the others posted: the additional bit of information you need is https://www.mathworks.com/help/matlab/math/parameterizing-functions.html

[Ekin Ceyda Cetin]

I was using griddedInterpolant to make the look up table. But I cannot use the griddedInterpolant when I calculate the objective function, because my input (value to look up) is a variable.

The problem is, during optimization when one parameter changes it should change a value in the objective function according to look up table.

Is there anything I can use besides griddedInterpolant? Or how to include this table in the optimization?

syms param1 param2
param=[param1,param2];

X.sort = sort(peakdata,'ascend');
X.r = size(X.sort);
X.rank = (1:X.r)';
CdfY = zeros(X.r(1),1);
CdfY = X.rank./(X.r(1)+1);
CdfX=X.sort;

a1=0;
a2=0;
a3=0;
for l=1:X.r(1)
a1=a1+log(CdfX(l)-param2);
a2=a2+(log((CdfX(l)-param2)*((CdfX(l)-param2)^param1)));
a3=a3+((CdfX(l)-param2)^param1);
end

b1=0;
b2=0;
b3=0;
for l=1:X.r(1)
b1=b1+((1/(CdfX(l)-param2)));
b2=b2+((CdfX(l)-param2)^param1);
b3=b3+((CdfX(l)-param2)^(param1-1));
end

xdata=linspace(0.2,1.2,1000);
vdata=z3;
lookup = griddedInterpolant(xdata,vdata);

%f3 is the objective function
f3=@(param) (z2/param1 +(1/X.r(1)*(a1))-(a2)/(a3))^2+(1/X.r(1)*(b1)*(b2)/(b3)-(lookup(param1)))^2;

low_b =[0,0];
up_b=[inf,peakdata(1)]
param0=[0.1,0.1];
options = optimoptions('fmincon', ...
    'Algorithm','interior-point', ...   
    'Display','iter-detailed');

[xfinal,fval,exitflag,output] = fmincon(f3, param0, [], [], [], [], low_b, up_b,[], options);

[Walter Roberson]

You are attempting to do interpolation at a symbolic location param1 and you are attempting to use an objective function that is defined symbolically, in terms of the symbolic variables param1 and param2 and the expressions a1, a2, a3, b1, b2, b3 that are defined in terms of those symbolic variables.

Try something more like,

syms param1 param2
param=[param1,param2];

X.sort = sort(peakdata,'ascend');
X.r = size(X.sort);
X.rank = (1:X.r)';
CdfY = zeros(X.r(1),1);
CdfY = X.rank./(X.r(1)+1);
CdfX=X.sort;

a1=0;
a2=0;
a3=0;
for l=1:X.r(1)
a1=a1+log(CdfX(l)-param2);
a2=a2+(log((CdfX(l)-param2)*((CdfX(l)-param2)^param1)));
a3=a3+((CdfX(l)-param2)^param1);
end

b1=0;
b2=0;
b3=0;
for l=1:X.r(1)
b1=b1+((1/(CdfX(l)-param2)));
b2=b2+((CdfX(l)-param2)^param1);
b3=b3+((CdfX(l)-param2)^(param1-1));
end

xdata=linspace(0.2,1.2,1000);
vdata=z3;
lookup = griddedInterpolant(xdata,vdata);

A1 = matlabFunction(a1, 'vars', {param(:)});
A2 = matlabFunction(a2, 'vars', {param(:)});
A3 = matlabFunction(a3, 'vars', {param(:)});
B1 = matlabFunction(B1, 'vars', {param(:)});
B2 = matlabFunction(B2, 'vars', {param(:)});
B3 = matlabFunction(B3, 'vars', {param(:)});

%f3 is the objective function
f3=@(param) (z2/param(1) +(1/X.r(1)*(A1(param)))-(A2(param))/(A3(param)))^2+(1/X.r(1)*(B1(param))*(B2(param))/(B3(param))-(lookup(param(1))))^2;

low_b =[0,0];
up_b=[inf,peakdata(1)]
param0=[0.1,0.1];
options = optimoptions('fmincon', ...
    'Algorithm','interior-point', ...   
    'Display','iter-detailed');

[xfinal,fval,exitflag,output] = fmincon(f3, param0, [], [], [], [], low_b, up_b,[], options);

[Ekin Ceyda Cetin]

This works perfectly, i didn't even know A1(param) was possible inside a function handle.

Another question, I suspect it is about the same thing. I need to use 'mu' and 'cdf(weibull)' in the below code. Therefore I put all the separate terms in the objective function at once using 'check_Yplot1'. It doesn't seem to work, the problem is parameter 'b' of the weibull distribution. Error:'The parameter B must be a positive finite numeric scalar.'

syms param1 param2;
param=[param1,param2];
X.r_check=size(peakdata);
check_XPlot1=peakdata;
mc5=linspace(1,X.r_check(1),X.r_check(1));

mu=(((1/(X.r_check(1)*w1))*(sum((check_XPlot1(mc5(1:X.r_check(1)))-param2).^param1)))^1/param1);
MU = matlabFunction(mu, 'vars', {param});
check_YPlot1 = @(param)cdf(makedist('weibull','a',MU(param),'b',param1),peakdata-param2);
f3=@(param)(((sum((log(1./(1-check_YPlot1(mc5(1:X.r_check(1))))).*log(log(1./(1-check_YPlot1(mc5(1:X.r_check(1)))))))))/(sum(log(1./(1-check_YPlot1(mc5(1:X.r_check(1)))))))-((1/X.r_check(1))*(sum(log(log(1./(1-check_YPlot1(mc5(1:X.r_check(1))))))))))/param1 +(1/X.r_check(1)*(sum(log(check_XPlot1(mc5(1:X.r_check(1)))-param2))))-(sum(log((check_XPlot1(mc5(1:X.r_check(1)))-param2).*((check_XPlot1(mc5(1:X.r_check(1)))-param2).^param1))))/(sum((check_XPlot1(mc5(1:X.r_check(1)))-param2).^param1)))^2+(1/X.r_check(1)*(sum((1./(check_XPlot1(mc5(1:X.r_check(1)))-param2))))*(sum((check_XPlot1(mc5(1:X.r_check(1)))-param2).^param1))/(sum((check_XPlot1(mc5(1:X.r_check(1)))-param2).^(param1-1)))-(((1/X.r_check(1))*(sum(log(1./(1-check_YPlot1(mc5(1:X.r_check(1))))))))*(sum((log(1./(1-check_YPlot1(mc5(1:X.r_check(1)))))).^(-1/param1)))/(sum((log(1./(1-check_YPlot1(mc5(1:X.r_check(1)))))).^((param1-1)/param1)))))^2;

low_b =[0,0];
up_b=[inf,peakdata(1)];
param0=[0.2,0.2];
options = optimoptions('fmincon', ...
    'Algorithm','interior-point', ...   
    'Display','iter-detailed');

[xfinal,fval,exitflag,output] = fmincon(f3, param0, [], [], [], [], low_b, up_b,[], options);

[Walter Roberson]

I notice that when you invoked matlabFunction that you used {param} as the vars. That creates a function handle that expects two separate inputs. However, fmincon passes the parameters as a single vector variable with multiple rows. The (:) that I used in defining the matlabFunction is important for telling matlabFunction to expect that the different inputs fall along the rows rather than the columns.

[Ekin Ceyda Cetin]

If I use it with (:), it says 'Index exceeds matrix dimensions'. If i don't use (:) problem with 'b' of weibull continues. Can it be related that I didn't supply a gradient function?

peakdata = load('experiment.50.p18.txt');
syms param1 param2;
param=[param1,param2];

mc5=linspace(1,X.r_check(1),X.r_check(1));

X.r_check=size(peakdata);
check_XPlot1=peakdata;
X.rank_check = (1:X.r_check)';
check_YPlot1 = zeros(X.r_check(1),1);
check_YPlot1 = X.rank_check./(X.r_check(1)+1); 

mu=((1/(X.r_check(1)*((1/X.r_check(1))*(sum(log(1./(1-check_YPlot1(mc5(1:X.r_check(1)))))))))*(sum((check_XPlot1(mc5(1:X.r_check(1)))-param2).^param1)))^1/param1);
MU = matlabFunction(mu, 'vars', {param(:)});
check_YPlot1 = @(param)cdf(makedist('weibull','a',MU(param),'b',param1),peakdata-param2);
f3=@(param)(((sum((log(1./(1-check_YPlot1(mc5(1:X.r_check(1))))).*log(log(1./(1-check_YPlot1(mc5(1:X.r_check(1)))))))))/(sum(log(1./(1-check_YPlot1(mc5(1:X.r_check(1)))))))-((1/X.r_check(1))*(sum(log(log(1./(1-check_YPlot1(mc5(1:X.r_check(1))))))))))/param1 +(1/X.r_check(1)*(sum(log(check_XPlot1(mc5(1:X.r_check(1)))-param2))))-(sum(log((check_XPlot1(mc5(1:X.r_check(1)))-param2).*((check_XPlot1(mc5(1:X.r_check(1)))-param2).^param1))))/(sum((check_XPlot1(mc5(1:X.r_check(1)))-param2).^param1)))^2+(1/X.r_check(1)*(sum((1./(check_XPlot1(mc5(1:X.r_check(1)))-param2))))*(sum((check_XPlot1(mc5(1:X.r_check(1)))-param2).^param1))/(sum((check_XPlot1(mc5(1:X.r_check(1)))-param2).^(param1-1)))-(((1/X.r_check(1))*(sum(log(1./(1-check_YPlot1(mc5(1:X.r_check(1))))))))*(sum((log(1./(1-check_YPlot1(mc5(1:X.r_check(1)))))).^(-1/param1)))/(sum((log(1./(1-check_YPlot1(mc5(1:X.r_check(1)))))).^((param1-1)/param1)))))^2;

low_b =[0,0];
up_b=[inf,peakdata(1)];
param0=[0.2,0.2];
options = optimoptions('fmincon', ...
    'Algorithm','interior-point', ...   
    'Display','iter-detailed');   
[xfinal,fval,exitflag,output] = fmincon(f3, param0, [], [], [], [], low_b, up_b,[], options);

[Walter Roberson]

You have

check_YPlot1 = @(param)cdf(makedist('weibull','a',MU(param),'b',param1),peakdata-param2);

Here you have bound the word "param" as a local variable for the purpose of anonymous function. Inside the @(param) body, "param" no longer means [param1,param2]. So there is no connection between the @(param) and the param1 and param2 in the body of the anonymous function. If you are expecting a vector of two parameters then you should index them:

check_YPlot1 = @(param)cdf(makedist('weibull','a',MU(param),'b',param(1)),peakdata-param(2));

This would create an anonymous function that expected a vector of two values as input and would call cdf with expressions based on the parameters, getting back a vector of results. That sounds like something you might plausibly want to do.

Now look for a moment at mc5

X.r_check=size(peakdata);
mc5=linspace(1,X.r_check(1),X.r_check(1));

That definition of mc5 is the same as

mc5 = 1 : X.r_check(1);

and establishes mc5 as an array of row indices, of length X.r_check(1). So mc5(1:X.r_check(1)) is going to be the same as mc5(1:end) which is going to be just 1 : X.r_check(1) -- consecutive integers. And there are going to be as many of them as your peakdata is tall.

Now look at your

f3=@(param)(((sum((log(1./(1-check_YPlot1(mc5(1:X.r_check(1))))).*log(log(1./(1-check_YPlot1(mc5(1:X.r_check(1)))))))))/(sum(log(1./(1-check_YPlot1(mc5(1:X.r_check(1)))))))-((1/X.r_check(1))*(sum(log(log(1./(1-check_YPlot1(mc5(1:X.r_check(1))))))))))/param1 +(1/X.r_check(1)*(sum(log(check_XPlot1(mc5(1:X.r_check(1)))-param2))))-(sum(log((check_XPlot1(mc5(1:X.r_check(1)))-param2).*((check_XPlot1(mc5(1:X.r_check(1)))-param2).^param1))))/(sum((check_XPlot1(mc5(1:X.r_check(1)))-param2).^param1)))^2+(1/X.r_check(1)*(sum((1./(check_XPlot1(mc5(1:X.r_check(1)))-param2))))*(sum((check_XPlot1(mc5(1:X.r_check(1)))-param2).^param1))/(sum((check_XPlot1(mc5(1:X.r_check(1)))-param2).^(param1-1)))-(((1/X.r_check(1))*(sum(log(1./(1-check_YPlot1(mc5(1:X.r_check(1))))))))*(sum((log(1./(1-check_YPlot1(mc5(1:X.r_check(1)))))).^(-1/param1)))/(sum((log(1./(1-check_YPlot1(mc5(1:X.r_check(1)))))).^((param1-1)/param1)))))^2;

You invoke that check_YPlot1 function multiple times, passing in to it mc5(1:X.r_check(1)). We just established above that mc5 expression is the same as 1:X.r_check(1) itself, the vector of integers. So you are calling check_YPlot1 with 1:X.r_check(1) as the vector of arguments, which will not typically be of length 2 as your peakdata will usually have more than 2 rows.

You are passing that vector of integers into check_YPlot1 that expects a vector of exactly two parameters, corresponding to param1 and param2. If that is deliberate, then you are effectively calling check_YPlot1([1,2]) a whole bunch of times. I can't see a good reason for that.

AGGGH

You have

check_YPlot1 = zeros(X.r_check(1),1);
check_YPlot1 = X.rank_check./(X.r_check(1)+1); 

mu=((1/(X.r_check(1)*((1/X.r_check(1))*(sum(log(1./(1-check_YPlot1(mc5(1:X.r_check(1)))))))))*(sum((check_XPlot1(mc5(1:X.r_check(1)))-param2).^param1)))^1/param1);

so you first define check_YPlot1 as a vector, then you overwrite it with what looks likely to be another vector, then you define mu in terms of that vector and subscripts, and then you define a function check_YPlot1. And then you define f3 in terms of a whole bunch of calls to that function, with arguments that look exactly like what you used when you were indexing a vector! This is getting really confusing!

I do see now that your f3 = @(param) uses param1 and param2, and perhaps you were thinking of those when you wrote your function check_YPlot1 ???

You need to define more clearly:

• what parameters are to be passed in to the function check_YPlot1 ?
• How does mc5(1:X.r_check(1)) get you those parameters?
• Are you trying to define your mu and check_YPlot1 recursively ??

Have another look at your f3 function. You have a bunch of calls to

check_YPlot1(mc5(1:X.r_check(1))

For the purpose of analysis, assign that value to a variable such as YP1 everywhere in a copy of the f3 definition.

YP1 = check_YPlot1(mc5(1:X.r_check(1));

and replace in the code. You will then see that you have a bunch of (1-YP1) expressions and that you have no other check_YPlot1 left. So simplify again, say calling it OneYP

OneYP = 1 - YP1;

You will now find that all OneYP occur in the form log(1./OneYP) so you can reduce that again,

logrOneYP = log(1./OneYP);

Your code is now much shorter and you finally do not always use that variable the same way. So, now back-substitute and turn it into a function:

Logr = @(param) log(1./(1-check_YPlot1(param)));

and write calls to that function in place of the longer expression in f3.

Now you still have references to check_XPlot1(mc5(1:X.r_check(1)) in the function and you could simplify those. The earlier analysis above suggests that can be replaced by just check_XPlot1 (the entire vector).

I will leave off there due to confusion about what you are really trying to calculate.

[Ekin Ceyda Cetin]

Thank you for your patience. This was enlightening information.

First of all, check_YPlot1(mc5(1:X.r_check(1)) should be just check_YPlot1(param),stupid of me. In the beginning, before overwriting it with function was to calculate mu for beginning. Now I got rid of that part where mu needs check_YPlot1 to be calculated. In the current code, I changed check_YPlot1, since I need to put my data in the function, it became YPlot(param,peakdata). Is there any other way to pass peakdata in there?

I made some simplifications in the f3 as you suggested with c1,c2,c3,d2,d3,a1,a2,a3,b1,b2,b3. However I cannot get it to work.

I don't know if it is ok to pass peakdata in the objective function as @(param,peakdata). Peakdata is a fixed row vector.

What I need to do is to let 'mu' change according to param1 and param2. Then use mu to get cdf vector from Yplot to calculate the objective function again to get next param1 and param2 iteratively. Is this possible?

peakdata = load('experiment.20.p19.txt');
syms param1 param2;
param=[param1,param2];
Xr=size(peakdata);
X.rank = (1:Xr(1))';
CdfY = zeros(Xr(1),1);
CdfY = X.rank./(Xr(1)+1);
d1=sum(log(1./(1-CdfY(1:Xr(1)))));
w1=(1/Xr(1))*(d1);
mu=(((1/(Xr(1)*w1))*(sum((peakdata(1:Xr(1))-param(2)).^param(1)))).^1/param(1));
MU = matlabFunction(mu, 'vars', {param(:)});
YPlot = @(param,peakdata)cdf(makedist('weibull','a',MU(param),'b',param(1)),peakdata-param(2));
c1= sum(log(1./(1-YPlot(param,peakdata))).*log(log(1./(1-YPlot(param,peakdata)))));
c2= sum(log(1./(1-YPlot(param,peakdata))));
c3= sum(log(log(1./(1-YPlot(param,peakdata)))));
%w2=(c1)/(c2)-((1/Xr(1))*(c3));
d2= sum(log(1./(1-YPlot(param,peakdata)))).^(-1/param(1));
d3= sum(log(1./(1-YPlot(param,peakdata)))).^((param(1)-1)/param(1));
%w3= w1*(d2)/(d3);
a1=sum(log(peakdata(1:Xr(1))-param(2)));
a2=sum(log((peakdata(1:Xr(1))-param(2)).*((peakdata(1:Xr(1))-param(2)).^param(1))));
a3=sum((peakdata(1:Xr(1))-param(2)).^param(1));
%f1=(w2/param1 +(1/X.r(1)*(A1(param)))-(A2(param))/(A3(param)))^2;
b1=sum((1./(peakdata(1:Xr(1))-param(2))));
b2=sum((peakdata(1:Xr(1))-param(2)).^param(1));
b3=sum((peakdata(1:Xr(1))-param(2)).^(param(1)-1));
%f2=(1/Xr(1)*(B1(param))*(B2(param))/(B3)-(w3))^2;
f3=@(param,peakdata)(((c1)/(c2)-((1/Xr(1))*(c3)))/param(1) +(1/Xr(1)*(a1))-(a2)/(a3))^2+(1/Xr(1)*(b1)*(b2)/(b3)-(w1*(d2)/(d3)))^2
%f3=@(param,peakdata)((w2)/param(1) +(1/Xr(1)*(a1))-(a2)/(a3))^2+(1/Xr(1)*(b1)*(b2)/(b3)-(w3))^2
%f3=f1+f2
low_b =[0,0];
up_b=[inf,peakdata(1)];
param0=[0.2,2];
options = optimoptions('fmincon', ...
    'Algorithm','interior-point', ... 
    'MaxFunctionEvaluations',100000,...
    'MaxIterations',80000,...
    'Display','iter-detailed');    
[xfinal,fval,exitflag,output] = fmincon(f3, param0, [], [], [], [], low_b, up_b,[], options);

Answer 2

[the cyclist]

The File Exchange has this code that says it functions like vlookup in Excel. You could try that.

[Ekin Ceyda Cetin]

Thank you for your answer. I need to you use interpolated values in the table as well so this code is not what I exactly need.

Answer 3

[Image Analyst]

Sure. If you have a finite set of input values but a ton of numbers to compute it on, then it makes sense to compute it just once for every possible input value and build a look up table. Then just apply the value from lookup table to every value. It will be much faster.

For example if you have some complicated formula

outputGrayLevel = SomeComplicatedFormula(inputGrayLevel);

and you need to apply that 20 million times for a 20 megapixel image, then it makes sense to just do it 256 times for every possible gray level (of an 8 but integer image) and then use

outputImage = intlut(inputImage, lut);

Using the built-in intlut() function will be much, much faster than computing that complicated formula 20 million times.