Why fmincon fails with this problem when the dimension increases?

Hi everyone,

I have a problem running the following code. I am trying to maximize a function called "sum rate" under nonlinear constraints given by the function "constraints" and when I run the main code using fmincon, most of the times it fails: either stops prematurely because of exceeding the function evaluation limit or the maximum iteration number, which I increased to 1e4 and 1e4 respectively) or it says that the function returned inf or -inf value though I am sure that it can't reach this value.

It is worth mentioning that the number of failed trials increase as I increase the dimensions of the problem (i.e. increase variable "M"). It works fair enough for M=3 but fails most of the time for M=5.

The problem is not convex by the way.

Here is the main code:

SNR = 10:10:100; %in dB
p = 10.^(SNR/10);
n = 2;
M = 2*n+1;
no_of_trials = 10;
options_fmin = optimset('Algorithm','active-set','FunValCheck','on','MaxFunEvals',1e4,'MaxIter',1e4);

for itrials=1:no_of_trials

    H11 = (randn(M,M)+sqrt(-1)*randn(M,M))/sqrt(2);   
    H12 = (randn(M,M)+sqrt(-1)*randn(M,M))/sqrt(2);
    H13 = (randn(M,M)+sqrt(-1)*randn(M,M))/sqrt(2);
    H21 = (randn(M,M)+sqrt(-1)*randn(M,M))/sqrt(2);
    H22 = (randn(M,M)+sqrt(-1)*randn(M,M))/sqrt(2);
    H23 = (randn(M,M)+sqrt(-1)*randn(M,M))/sqrt(2);
    H31 = (randn(M,M)+sqrt(-1)*randn(M,M))/sqrt(2);
    H32 = (randn(M,M)+sqrt(-1)*randn(M,M))/sqrt(2);
    H33 = (randn(M,M)+sqrt(-1)*randn(M,M))/sqrt(2);
    T = inv(H21)*H23*inv(H13)*H12*inv(H32)*H31;
    cond_no(itrials) = cond(T);
    [V D] = eig(T);
    d=diag(D);
    Gamma1(:,1) = ones(M,1);
    for i=1:n
        Gamma1(:,i+1) = d.^i;
        Gamma3(:,i) = d.^i;
    end
    Gamma2(:,1) = ones(M,1);
    for i=1:n-1
        Gamma2(:,i+1) = d.^i;
    end

        F11 = [];
        F12 = [];
        F13 = [];
        F21 = [];
        F22 = [];
        F23 = [];
        F31 = [];
        F32 = [];
        F33 = [];
        for i=1:M
            F11 = [F11; diag(H11(i,:)*V)*Gamma1];
            F12 = [F12; diag(H12(i,:)*inv(H32)*H31*V)*Gamma2];
            F13 = [F13; diag(H13(i,:)*inv(H23)*H21*V)*Gamma3];
            F21 = [F21; diag(H21(i,:)*V)*Gamma1];
            F22 = [F22; diag(H22(i,:)*inv(H32)*H31*V)*Gamma2];
            F23 = [F23; diag(H21(i,:)*V)*Gamma3];
            F31 = [F31; diag(H31(i,:)*V)*Gamma1];
            F32 = [F32; diag(H31(i,:)*V)*Gamma2];
            F33 = [F33; diag(H33(i,:)*inv(H23)*H21*V)*Gamma3];
        end

        M1n = [F11 sqrt(2)*F12];
        M1d = sqrt(2)*F12;
        P2 = [1 zeros(1,n); zeros(n,1) sqrt(2)*eye(n)];
        M2n = [F22 F21*P2];
        M2d = F21*P2;
        P3 = [sqrt(2)*eye(n) zeros(n,1); zeros(1,n) 1];
          M3n = [F33 F31*P3];
          M3d = F31*P3;

          w0_cmpx = randn(2*M,1);
          [w_fmin fval_fmin exitflag_fmin(itrials,isnr)]=fmincon(@(w_fmin)sum_rate_func_wcmpx_MxM_111312(w_fmin,V,M,n,p(isnr),M1n,M2n,M3n,M1d,M2d,M3d),w0_cmpx,[],[],[],[],[],[],@(w_fmin)constraint_func_wcmpx_MxM_111312(w_fmin,d,M,n,V,H21,H23,H31,H32,3*M),options_fmin);

end
end

Here is the sum rate function (the objective function):

function f=sum_rate_func_wcmpx_MxM_111312(w,V,M,n,p,M1n,M2n,M3n,M1d,M2d,M3d)

for i=1:M
    w_cmpx(i,1)=w(2*(i-1)+1) + sqrt(-1)*w(2*i);
end
w_tilda = inv(V)*w_cmpx;
W_block(1,:) = [w_tilda.' zeros(1,M^2-M)];
for i=2:M
    W_block(i,:) = circshift(W_block(i-1,:),[0,M]);
end
W_tilda_block = W_block'*W_block;

rate_1 = log2(det(eye(M)+ p*M1n'*W_tilda_block*M1n)) - log2(det(eye(n) + p*M1d'*W_tilda_block*M1d));
rate_2 = log2(det(eye(M)+ p*M2n'*W_tilda_block*M2n)) - log2(det(eye(n+1) + p*M2d'*W_tilda_block*M2d));
rate_3 = log2(det(eye(M)+ p*M3n'*W_tilda_block*M3n)) - log2(det(eye(n+1) + p*M3d'*W_tilda_block*M3d));

f = -real(rate_1+rate_2+rate_3);

Here is the nonlinear constraint function:

function [c,ceq] =   constraint_func_wcmpx_MxM_111312(w,d,M,n,V,H21,H23,H31,H32,con)

    for i=1:M
        w_cmpx(i,1)=w(2*(i-1)+1) + sqrt(-1)*w(2*i);
    end
    w_tilda = inv(V)*w_cmpx;
    Wd = diag(w_tilda);

    Gamma1(:,1) = ones(M,1);
    for i=1:n
        Gamma1(:,i+1) = d.^i;
        Gamma3(:,i) = d.^i;
    end
    Gamma2(:,1) = ones(M,1);
    for i=1:n-1
        Gamma2(:,i+1) = d.^i;
    end

    v1 = V*Wd*Gamma1;
    v2 = inv(H32)*H31*V*Wd*Gamma2;
    v3 = inv(H23)*H21*V*Wd*Gamma3;

    c = [];
    ceq = [real(trace(v1'*v1) + trace(v2'*v2) + trace(v3'*v3)) - (con)];