How to actually maximize a function with many variables

Hello. I have these two functions:

lik = @(x) -(log(x(1,1))*mb2(1,1)+log(x(1,2))*mb2(1,2)+log(x(1,3))*mb2(1,3)+log(x(1,4))*mb2(1,4)...
            +log(x(2,1))*mb2(2,1)+log(x(2,2))*mb2(2,2)+log(x(2,3))*mb2(2,3)+log(x(2,4))*mb2(2,4)...
            +log(x(3,1))*mb2(3,1)+log(x(3,2))*mb2(3,2)+log(x(3,3))*mb2(3,3)+log(x(3,4))*mb2(3,4)...
            +log(x(4,1))*mb2(4,1)+log(x(4,2))*mb2(4,2)+log(x(4,3))*mb2(4,3)+log(x(4,4))*mb2(4,4));

lik2 = @(y) -(log(y(1,1))*mb(1,1)+log(y(1,2))*mb(1,2)+log(y(1,3))*mb(1,3)+log(y(1,4))*mb(1,4)...
             +log(y(2,1))*mb(2,1)+log(y(2,2))*mb(2,2)+log(y(2,3))*mb(2,3)+log(y(2,4))*mb(2,4)...
             +log(y(3,1))*mb(3,1)+log(y(3,2))*mb(3,2)+log(y(3,3))*mb(3,3)+log(y(3,4))*mb(3,4)...
             +log(y(4,1))*mb(4,1)+log(y(4,2))*mb(4,2)+log(y(4,3))*mb(4,3)+log(y(4,4))*mb(4,4)...
             +log(y(1,5))*mb(1,5)+log(y(1,6))*mb(1,6)+log(y(1,7))*mb(1,7)+log(y(1,8))*mb(1,8)...
             +log(y(2,5))*mb(2,5)+log(y(2,6))*mb(2,6)+log(y(2,7))*mb(2,7)+log(y(2,8))*mb(2,8)...
             +log(y(3,5))*mb(3,5)+log(y(3,6))*mb(3,6)+log(y(3,7))*mb(3,7)+log(y(3,8))*mb(3,8)...
             +log(y(4,5))*mb(4,5)+log(y(4,6))*mb(4,6)+log(y(4,7))*mb(4,7)+log(y(4,8))*mb(4,8));

Where mb and mb2 are matrices of positive integers. The solution to this, x and y, are matrices of numbers between 0 and 1, but it is very important to me to be as accurate as possible (the logs are there because of that), and I am only getting 'Local minimum possible' when using fmincon, instead of 'Local minimum found'. I change x0 and y0 and sometimes get a better solution, but it still does not show 'Local minimum found'. Here are the constraints and the fmincom I am using:

    % Maximizing with constraints:

    x0 = mb2*(1/100);
    Aeq = [1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0; ... % Make each line of x sum to one.
           0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0; ...
           0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0; ...
           0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1];
    beq = ones(4,1);
    lb = zeros(4,4);
    ub = ones(4,4);
    nonlcon = @cons_am;
    options = optimoptions(@fmincon,'steptolerance',1.0000e-13,'constraintTolerance',1.0000e-09,'optimalitytolerance',1.0e-16,'MaxFunctionEvaluations',100000);
    [x, xval, exitflag] = fmincon(lik,x0,[],[],Aeq,beq,lb,ub,nonlcon,options);

    y0 = mb*(1/100);
    Aeq2 = [1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;...
            0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;...
            0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;...
            0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;...
            0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0;...
            0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0;...
            0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0;...
            0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1];

    beq2 = ones(8,1);
    lb2 = zeros(4,8);
    ub2 = ones(4,8);
    nonlcon2 = @cons_am2;
    options2 = optimoptions(@fmincon,'steptolerance',1.0000e-13,'constraintTolerance',1.0000e-09,'optimalitytolerance',1.0e-16,'MaxFunctionEvaluations',100000);
    [y, yval, exitflag2] = fmincon(lik2,y0,[],[],Aeq2,beq2,lb2,ub2,nonlcon2,options2);

function [c, ceq] = cons_am2(y)
mu = [0.032051282051282  0.0000;0.33974358974359  0.140127388535032;0.455128205128205  0.407643312101911;0.173076923076923  0.452229299363057];
c(1)=3-(0.5873*(mu(1,1)*y(1,1)+mu(2,1)*y(2,1)+mu(3,1)*y(3,1)+mu(4,1)*y(4,1))/(mu(1,2)*y(1,5)+mu(2,2)*y(2,5)+mu(3,2)*y(3,5)+mu(4,2)*y(4,5)));
c(2)=  (0.5873*(mu(1,1)*y(1,2)+mu(2,1)*y(2,2)+mu(3,1)*y(3,2)+mu(4,1)*y(4,2))/(mu(1,2)*y(1,6)+mu(2,2)*y(2,6)+mu(3,2)*y(3,6)+mu(4,2)*y(4,6)))-3;
c(3)=1-(0.5873*(mu(1,1)*y(1,2)+mu(2,1)*y(2,2)+mu(3,1)*y(3,2)+mu(4,1)*y(4,2))/(mu(1,2)*y(1,6)+mu(2,2)*y(2,6)+mu(3,2)*y(3,6)+mu(4,2)*y(4,6)));
c(4)=  (0.5873*(mu(1,1)*y(1,3)+mu(2,1)*y(2,3)+mu(3,1)*y(3,3)+mu(4,1)*y(4,3))/(mu(1,2)*y(1,7)+mu(2,2)*y(2,7)+mu(3,2)*y(3,7)+mu(4,2)*y(4,7)))-1;
c(5)=(1/3)-(0.5873*(mu(1,1)*y(1,3)+mu(2,1)*y(2,3)+mu(3,1)*y(3,3)+mu(4,1)*y(4,3))/(mu(1,2)*y(1,7)+mu(2,2)*y(2,7)+mu(3,2)*y(3,7)+mu(4,2)*y(4,7)));
c(6)=  (0.5873*(mu(1,1)*y(1,4)+mu(2,1)*y(2,4)+mu(3,1)*y(3,4)+mu(4,1)*y(4,4))/(mu(1,2)*y(1,8)+mu(2,2)*y(2,8)+mu(3,2)*y(3,8)+mu(4,2)*y(4,8)))-(1/3);
 ceq = [];
end

function [c, ceq] = cons_am(x)
mu = [0.032051282051282  0.0000;0.33974358974359  0.140127388535032;0.455128205128205  0.407643312101911;0.173076923076923  0.452229299363057];
c(1)=3-(0.5873*(mu(1,1)*x(1,1)+mu(2,1)*x(2,1)+mu(3,1)*x(3,1)+mu(4,1)*x(4,1))/(mu(1,2)*x(1,1)+mu(2,2)*x(2,1)+mu(3,2)*x(3,1)+mu(4,2)*x(4,1)));
c(2)=  (0.5873*(mu(1,1)*x(1,2)+mu(2,1)*x(2,2)+mu(3,1)*x(3,2)+mu(4,1)*x(4,2))/(mu(1,2)*x(1,2)+mu(2,2)*x(2,2)+mu(3,2)*x(3,2)+mu(4,2)*x(4,2)))-3;
c(3)=1-(0.5873*(mu(1,1)*x(1,2)+mu(2,1)*x(2,2)+mu(3,1)*x(3,2)+mu(4,1)*x(4,2))/(mu(1,2)*x(1,2)+mu(2,2)*x(2,2)+mu(3,2)*x(3,2)+mu(4,2)*x(4,2)));
c(4)=  (0.5873*(mu(1,1)*x(1,3)+mu(2,1)*x(2,3)+mu(3,1)*x(3,3)+mu(4,1)*x(4,3))/(mu(1,2)*x(1,3)+mu(2,2)*x(2,3)+mu(3,2)*x(3,3)+mu(4,2)*x(4,3)))-1;
c(5)=(1/3)-(0.5873*(mu(1,1)*x(1,3)+mu(2,1)*x(2,3)+mu(3,1)*x(3,3)+mu(4,1)*x(4,3))/(mu(1,2)*x(1,3)+mu(2,2)*x(2,3)+mu(3,2)*x(3,3)+mu(4,2)*x(4,3)));
c(6)=  (0.5873*(mu(1,1)*x(1,4)+mu(2,1)*x(2,4)+mu(3,1)*x(3,4)+mu(4,1)*x(4,4))/(mu(1,2)*x(1,4)+mu(2,2)*x(2,4)+mu(3,2)*x(3,4)+mu(4,2)*x(4,4)))-(1/3);
ceq = [];
end

If anyone has any suggestions I would be very grateful. I am buffled for not being able to get it to work properly. In case some context works: I have two samples from an experiment, I am trying to choose between two hypothesis by running a Neyman-Pearson test. The functions to maximize are the maximum likelihoods, and the constraints are there to restrict the likelihood to some intervals of my interest. Thanks again.

[OCDER]

Whoa, what are we looking at o_O? Your code is more like an ASCII art like this:

("`-''-/").___..--''"`-._ 
 `6_ 6  )   `-.  (     ).`-.__.`) 
 (_Y_.)'  ._   )  `._ `. ``-..-' 
   _..`--'_..-_/  /--'_.'
  ((((.-''  ((((.'  (((.-' 

 Only, this ASCII art is simpler.

Isn't the following :

lik2 = @(y) -(log(y(1,1))*mb(1,1)+log(y(1,2))*mb(1,2)+log(y(1,3))*mb(1,3)+log(y(1,4))*mb(1,4)...
             +log(y(2,1))*mb(2,1)+log(y(2,2))*mb(2,2)+log(y(2,3))*mb(2,3)+log(y(2,4))*mb(2,4)...
             +log(y(3,1))*mb(3,1)+log(y(3,2))*mb(3,2)+log(y(3,3))*mb(3,3)+log(y(3,4))*mb(3,4)...
             +log(y(4,1))*mb(4,1)+log(y(4,2))*mb(4,2)+log(y(4,3))*mb(4,3)+log(y(4,4))*mb(4,4)...
             +log(y(1,5))*mb(1,5)+log(y(1,6))*mb(1,6)+log(y(1,7))*mb(1,7)+log(y(1,8))*mb(1,8)...
             +log(y(2,5))*mb(2,5)+log(y(2,6))*mb(2,6)+log(y(2,7))*mb(2,7)+log(y(2,8))*mb(2,8)...
             +log(y(3,5))*mb(3,5)+log(y(3,6))*mb(3,6)+log(y(3,7))*mb(3,7)+log(y(3,8))*mb(3,8)...
             +log(y(4,5))*mb(4,5)+log(y(4,6))*mb(4,6)+log(y(4,7))*mb(4,7)+log(y(4,8))*mb(4,8));

the same as this??

lik2 = @(y) -sum(log(y(:).*mb(:)));

Before we can offer help for your main problem, perhaps ask the forum:

"How can I simplify my code? It's too complex..."

[Juan Gambetta]

Haha, well, thank you. I am a true beginner, and writing all that long ugly thing helps me see the problem more clearly. Have you got any ideas on my main issue?

[OCDER]

Hm... TBH, I can't even figure out what you're trying to do/solve because the code is too complex. Try to simplify, otherwise, change the question to be about simplifying your code first. Once that's fixed, then someone can help with the real problem.

EXAMPLE:

"How can I simplify my code? It's too complex..."

I'm trying to maximize my function by changing variables X, BUT, I've been told the code is too complex. Can someone help me reduce these complex codes? I'm new and my CTRL, C, and V keys can only take so much damage... Thanks!

lik = .... %any simpler way

like = .... %any simpler way?

function cons_am ... %any simpler way

function cons_am2 ... %any simpler way

[Juan Gambetta]

Basically, I want to minimize:

lik = @(x) -sum(x(:).*mb2(:))
lik2 = @(y) -sum(y(:).*mb(:));

Knowing both x (4x4 matrix) and y (4x8 matrix) are numbers between 0 and 1, and mb and mb2 are matrices of positive integers. The thing is, fmincon only returns 'Local minimum possible.' That is not enough, since I change the initial guess and the solution improves. How can I obtain the best solution there is?

I have seriously started to simplify my code in the meantime, but that is a question I have tried to answer for days.

[Stephen23]

"...writing all that long ugly thing helps me see the problem more clearly."

Learn to write simpler code, such as vectorized code that OCDER showed you. Complex code has more bugs and is harder to debug.

A good rule of thumb: any time you copy-and-paste code then you are doing something wrong.

[Juan Gambetta]

Thank you. I appreciate the advice.