MATLAB: Fmincon gives different optimal vector after adding equality constraint that describes output of fmincon before adding the constraint

fminconoptimizationOptimization Toolbox

Hello guys, I'm new to MathWorks but I think I have a decent experience with MATLAB.
So, I have a function called MaxDiversification which, in addition to other inputs, takes a structured array of inequality and equality constraints for my problem. When I run the code I get a vector, w, of weights for a portfolio, and the portfolio return, R, for that set vector of weights, i.e. R = Mu'*w where Mu is a vector of expected returns for each asset.
Now here's my problem: when I add the constraint Mu'*w = R to the structured array… where R is the output from before… and then run MaxDiversification, I get a different weights vector w. WHY?
I'm adding this constraint because I want to test the function before running it with other portfolio returns R and getting the corresponding optimal vectors w.
The code for the function is given below. I believe ENB and the functions MinimumTorsion and PrincipalComp don't affect the optimization of w so I didn't mention anything about them.
 function [w, ENB, MinRet] = MaxDiversification2(Mu, w_0, S, Constr, method)
switch method
    case 'MT'
        w = fmincon(@exp_term1, w_0, Constr.A, Constr.b, Constr.Aeq, Constr.beq);
        ENB = MinimumTorsion(w,S);
        MinRet = Mu'*w;
    case 'PCA'
        w = fmincon(@exp_term2, w_0, Constr.A, Constr.b, Constr.Aeq, Constr.beq);
        ENB = PrincipalComp(w, S);
        MinRet = Mu'*w;
end
    % exp_term1
    function y = exp_term1(w)
    Sigma = S;
    volatilities = diag(Sigma).^(1/2);
    corr = diag(1./volatilities)*Sigma*diag(1./volatilities);
    [eigvec,eigval] = eig(corr);
    gamma = sqrt(eigval);
    c = eigvec*gamma*eigvec';
    n = size(Sigma, 2);
    d = eye(n);
    pi_1 = eye(n);
    abserror_pi = norm(d,'fro');
    niter = 0;
    while abserror_pi > 10^-14
        niter = niter + 1;
        if niter == 8000
            break
        end
    pi_0 = pi_1;
    u = (d*(c*c')*d)^0.5;
    q = (u\d)*c;
    d = diag(diag(q*c));
    pi_1 = d*q;
    abserror_pi = norm(pi_1- pi_0,'fro')/n;
    end
    pi = pi_1;
    clear pi_1
    t_MT = diag(volatilities)*pi*(c\diag(1./volatilities));
    dist_MT = ((t_MT'\w).*(t_MT*Sigma*w))/(w'*Sigma*w);
    y = -exp(-dist_MT'*log(dist_MT));
    end
    % exp_term2
    function y = exp_term2(w)
    load S
    Sigma = S;
    [eigvec, eigval] = eig(Sigma);
    [~, ind] = sort(diag(eigval),'descend');
    eigvec = eigvec(:,ind);
    flip = eigvec(1,:) < 0;
    eigvec(:,flip) = -eigvec(:,flip);
    dist_PC = (eigvec'*w).*(eigvec'*Sigma*w)/(w'*Sigma*w);
    y = -exp(-dist_PC'*log(dist_PC));
    end
end   
  end

Best Answer

I can think of a a few possible reasons,
  1. Floating point arithmetic? You haven't said how big or small the differences are.
  2. One or both of the optimizations didn't actually complete. I notice there is not exitflag check in your code.
  3. Your solution space is non-unique. Did you compare the objective function value reached and whether constraints are satisfied for both cases?
  4. You have local minima and one or both optimizations got stuck at sub-optimal points.
Related Question